
QUASI-OPTIMUM DESIGN OF 
A SIX DEGREE OF FREEDOM 
MOVING BASE SIMULATOR 
CONTROL SYSTEM 

by Bernard Friedland, Chong-Kuan Ling, 
and Maurice F. Hutton 


Prepared by 
SINGER COMPANY 
Little Falls, N.J. 
for Ames Research Center 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION • WASHINGTON, D. C. • OCTOBER 1973 



1. Report No. 2. Government Accession No. 

NASA CR-2312 

3. Recipient's Catalog No. 

4. Title and Subtitle 

Quasi-Optimum Design of a Six Degree of Freedom Moving Base 
Simulator Control System 

5. Report Date 
October 1973 

6. Performing Organization Code 

7. Author(s) 

Bernard Friedland, Chong-Kuan Ling, and Maurice F. Hutton 

8. Performing Organization Report No. 

10. Work Unit No. ^ 

9. Performing Organization Name and Address 
Singer Company 
Kearfott Division 
Little Falls, New Jersey 

11. Contract or Grant No. 
NAS 2-3636 

13. Type of Report and Period Covered 
Contractor Report 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 

14. Sponsoring Agency Code 


15. Supplementary Notes 


16. Abstract 


The design of a "washout" control system for a moving-base simulator is treated by a quasi-optimum 
control technique developed by Firedland (NASA CR-527 and CR-1099) . An earlier version of the 
washout system, limited to longitudinal motion, was described in NASA 01-1613 and this work extends 
the design to a complete six -degree -of -freedom simulator. The broad objective of the design is to 
reproduce the sensed motion (angular velocity and specific force) as accurately as possible without 
causing the simulator (cab) excursions to exceed specified limits. A performance criterion is 
established that weights magnitude and direction errors in specific force and in angular velocity 
and attempts to maintain the excursion within set limits by penalizing excessive excursions. A 
FORTRAN routine for realizing the washout law was developed and typical time-histories using the 
washout routine were simulated for a range of parameters in the penal ty-and weighting -functions. 
These time histories and the listing of the routine are included in the report. Evaluation of 
this class of washout laws in actual pilot simulations was not undertaken, and such evaluation is 
recommended for future investigation. 


17. Key Words (Suggested by Author(s)) 

Simulator, Moving -base simulator, Optimal Control! 
Washout Control Law 


18. Distribution Statement 


UNCLASSIFIED-UNLIMITED 


19. Security Classif. (of this report) 
UNCLASSIFIED 

20. Security Classif. (of this page) 

UNCLASSIFIED 

21. No. of Pages 
149 . 

22. Price* 


- - - 


* For sale by the National Technical Information Service, Springfield, Virginia 22151 




















PREFACE. 


An earlier report [ 1 ] describes the application of a quasi-optimum 
control technique to the design of a control system for a three degree-of- 
freedom motion simulator. This technique, developed under this contract, 
NAS 2-3636 and an earlier contract, NAS 2-3648, was applied to design 
a complete six degree-of-freedom motion simulation in the investigation 
described herein. This report contains the analytical results and simulated 
time histories that would be obtained with various parameter settings for 
several types of missions. To facilitate experimental evaluation of the 
control law, a description and listing of the washout subroutine Is given 
in Appendix IV. 

The authors are grateful for the assistance provided by Mr. J. G. 
Douvillier and by Dr. E. C. Stewart, who served as Contract Technical 
Monitor. 
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1. INTRODUCTION 


In controlling the motion of an aircraft, the pilot uses both visual and kinesthetic 

a 

cues. Visual cues provide information about the position of the aircraft with respect 
to a suitable frame of reference; kinesthetic cues provide information about the motion , 
(velocity, acceleration) of the aircraft with respect to the reference frame. The 
pilot processes the data from the visual and kinesthetic sensors in a sophisticated (and 
little understood) manner in such a way that his action upon the aircraft controls , 
causes the aircraft to behave as he desires. 

Since the primary objective of the pilot is to control the position of the aircraft, 
it is reasonable to suppose that the pilot relies primarily upon visual cues; kinesthetic 
cues, although very useful, are not indispensible. Presumably, this is the reason that 
fixed-base (non-moving) simulators, which provide accurate visual cues, (but no kinesthetic 
cues) have been used with great success in a variety of aircraft and space applications. 

Owing to the absence of motion cues in fixed-base simulators, however, the ex- 
perience of the pilot in such simulators is not identical to what he would experience in 
an actual aircraft, and hence there are many instances in which the pilot's performance 
is not the same as it would be in an actual aircraft, and in which he may complain 

ft 

of the lack of fidelity of the simulation. It is generally believed that in the absence of 
motion cues, the pilot after an initial training period in the simulator, alters his method* 
of mental data processing in an attempt to maintain his performance at the level he would 
achieve in the actual aircraft. In effect, the pilot synthesizes or interpolates the in- 
formation he expects from his kinesthetic sensors. To accomplish this, however, requires 
greater mental effort; accordingly, the pilot becomes more rapidly fatigued, and he 
generally describes the task as being more difficult to perform. 

The recognized shortcomings of fixed-base simula tors hove led to the increas ing 
use of simulators which can move in response to the pilot's commands. The introduction 
of additional motion into simulators, however, has not alleviated all the problems of fixed- 
base simulators. There are even situations in which the pilot indicates preference for a 
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fixed-base simulation over a moving-base simulation. The major difficulty with 
moving-base simulators is that their motion is generally not identical to that of 
the aircraft being simulated, because the simulator is confined to a physical 
volume which is much smaller than the volume in which the aircraft is free to maneuver 
(The exception to this general situation is the so-called case of "one-to-one" motion, 
in which the pilot's task and aircraft are selected so that the aircraft could be maintained 
within a volume not larger than the volume available for motion of the simulator. An 
example of such a situation is a helicopter with a hovering task.) 

It is evident that fidelity of the motion cues can be increased by increasing 
the volume in which the simulator is free to maneuver. The important question in regard 
to design of motion simulators is thus how to make most effective use of a given 
maneuverability volume. 

To make effective use of the volume in which the simulator moves , some knowledge 
of the nature of the kinesthetic sensors of motion is required. It is generally believed 
that the principal sensors of motion in the human being are in the labyrinth structure of 
the ear, and comprise the semicircular canals and the otoliths. The former are believed 
to act in the manner of rate gyros to sense the angular velocity of the pilot's head and 

■ 

.the latter act like linear accelerometers to sense the specific force at the pilot's head. Motion 
of the body is also sensed as a result motion of the organs in the abdominal cavity and due to 
pressure on the body surface (the "seat of the pants"), but it is believed that the motion cues 
derived from other than the ear labyrinth are relatively insignificant. 

Various studies have been made to determine the characteristics of the labyrinthine 
sensors and it has been concluded that these are fairly complex nonlinear dynamic systems . 
Moreover, the characteristics of the physical sensors themselves, even if obtainable, 
would not be adequate to determine how the human being senses motion, because the 
signals from the sensors are processed by the sophisticated digital computer which is the 
brain. Sensor models have been constructed to account for the signal processing in the 
brain, but these are not very well established. Hence, for the purpose of this investigation, 
we have assumed that the motion cues of significance are simply angular velocity and 
specific force. Each are vector quantities and are refered to a set of reference axes fixed 
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in the pilot and moving with him. 

. ■*!_•,* • ' • .... 

The ideal motion simulator would produce in ti.e simulator cab the same angular 
velocity vector and specific force vector as would result in the aircraft in response 
to the pilot's control inputs. The difference between the angular velocity vector of 
the simulator and the angular velocity vector in' the aircraft is an error as is the dif- 
ference between the specific force vector in the cab and in the aircraft. If these- 
error vectors are both zero, the motion simulation is perfect ("one-to-one"). The 
objective of the simulator control system design is to keep these errors as small as 
possible. 

One of the difficulties in designing an optimum simulator control system is that there 
are two (2) vector-valued error components (angular velocity and specific force) or a 
total of 6 scalar error signals. It is thus necessary to devise a single scalar measure 
of error to account for the 6 components which may be present. Since the typical 
pilot is usually able to discriminate between types of motion simulations and can assess 
their relative performance, it is possible that such a measure of error exists. If it does 
exist, it may depend on the particular task, the pilot's experience, the nature of the 
visual cues, and many other factors. Hence determination of the "natural" measure 
of error, i.e. the measure of error which governs the pilot's subjective evaluation of 
simulator quality, is unrealistic. The only feasible approach is to use a "reasonable" 
measure of error’. The approach adopted in this study is to regard the total measure of 
error as the weighted sum of two scalars, one representing specific force, and one re- 
presenting angular velocity, i.e. 

E = M + k M 

3 co .... , 

i 

where M is the contribution to the total error due to errors in the specific force vector, 

P 

'M is the contribution due to errors in the angular velocity vector and k is a constant 

CO I — 

whichrdetermines the importance of angular velocity relative to specific force. The 
scalars and M are determined from the corresponding vectors on the basis of the 
following empirical considerations: 
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. The direction of the error is significant. For example, if the true 
specific force is 0.5 g and the simulated specific force is 1.5 g 
(i.e. the error magnitude is 1 g), this is less serious than when the 
simulated specific force is - 0.5 g (i.e.^the error is still 1 g in 
magnitude but the simulated specific force is in the opposite direction 
to the true specific force). 

. The specific force in the vertical direction when the aircraft is not 
accelerating vertically is 1 g (which is also the specific force ex- 
perienced by a person at rest on the ground). One is not normally 
aware of this vertical specific force when standing. In other words, 
it is reasonable to assume that the brain "biases-out" the normal 1 g 
component of vertical specific force, and hence that the vertical 
specific force sensation is the actual specific force less 1 g. 

Based on these considerations several alternative analytical expressions for the 
measure of error in specific force and angular velocity have been developed and are 
described in Section 2. 4- 

Using these error measures and the methods of optimum control theory, a set of 
control laws were devised. 
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2. ANALYSIS 


2.1 Description of Simulation Problem 

In the operation of a moving-base flight simulation, a pilot manipulates a set of flight 
controls and the simulator cab in which he is situated moves in a manner which tends to re- 
produce^ the sensation of motion which the pilot would experience if he were in the actual 
aircraft and he manipulated the flight controls in the same manner. Ideally, the motion of 
the cab should be identical to that of the aircraft. Under conditions which permit the cab 
motion to reproduce the aircraft motion perfectly ("one-to-one" simulation) the sensation of 
motion in the simulator is the same as the sensation of motion in the actual aircraft. Generally, 
however, simulators are used to simulate flight tasks in which one-to-one simulation is physically 
impossible because the cab is confined to remain inside a fairly small physical volume. The 
simulator control system thus must be designed, not to reproduce the motion of the aircraft 
being simulated, since this is impossible, but rather to provide a sensation of motion which is 
as close as possible to the sensation of motion in the actual aircraft, but without causing the 
cab to exceed the physical limits of its motion. 


The moving-base simulator of concern in this investigation is of the type currently in use 
at the Ames Research Center; the cab is built into a gimbal system which provides three 
rotational degrees of freedom and the whole gimbal assembly is in turn provided with three 
translational degrees of freedom by a system of three mutually perpendicular tracks. (See 
Figure 2-1). Each degree of freedom of the cab motion can be independently controlled by 
the cab arive system. A digital computer is used to simulate the dynamics of the aircraft 
which permits the dynamics to be modeled to any degree of accuracy desired. The inputs to 
this computer are the flight commands resulting from the pilot's manipulation of the flight 
controls in the cab. The computer outputs are the state variables (velocity, angular rotation, 
and etc. ) describing the motion of the aircraft which in turn are used by the simulator control 
to co mpute_the.signals-needed-to-drive-the cab~motidnT If is the design of the simulator control 
that is the subject of this study. Figure 2-2 shows a functional block diagram of the overall 
simulator system. 
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Figure 2-1 

The Ames All-Axis Motion Simulator 
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Figure 2-2 

Functional Control System Block Diag 




Since each degree of freedom of the cab is driven by an electromechanical system 
which is designed as a high-speed position servo, it is reasonable to assume that the position 
(i.e. gimbal angles for rotational degrees of freedom, translations along tracks for trans- 
lational degrees of freedom) of the cab is identical to the drive system position commands, 
so long as the latter do not exceed the motion limits. If the motion limits are exceeded, 
mechanical switches are operated to engage a system designed to arrest the cab motion safely. 
When the safety switches are engaged, the simulation is ended, the cab is returned to a 
neutral position, and, after various safety checks are made, the simulation is again initiated. 
To avoid engaging the safety switches, the drive system has an electrical system designed to 
anticipate the engagement of the mechanical safety system and to cutoff power before cab 
motion causes the latter to be engaged. Accordingly, from the operational viewpoint, the 
simulator control should be designed so that the cab drive command signals do not cause either 
the electrical or mechanical safety limits to.be exceeded. 

Because of the physical construction of the simulator, the position of the cab at any 
instant of time is approximately described by the three translations of the movable gimbal 
structure on the system of rails, and three gimbal angles. These quantities and their time 
derivatives (twelve in all) are a set of state variables natural to the motion of the cab. They 
are, however, not particularly well-suited for describing the sensible motion of the aircraft. 

For the latter, a more suitable set of state variables are the components of the vehicle linear 
and angular velocity vectors resolved along a set of axes fixed in the aircraft at the pilot's 
station (6 quantities) and 3 angles relating the position of these body axes to a set of reference 
axes. With regard to the sensation of motion, the position in space of the aircraft is not im- 
portant. For that matter, velocity, per se, is not important; it is acceleration that is re- 
sponsible for sensation of motion: linear translation at constant velocity does not contribute 

to the sensation of motion. This fact is particularly convenient with regard to simulation of 
the forward motion of the aircraft. Under quescent conditions (cruise, for example) the air- 
craft moves at constant forward velocity which is not sensible, except through visual reference. 
Only changes from this quiescent state are detectable; consequently, the constant forward 
component of velocity is not required in the simulation. As a consequence, the simulation is 
accomplished by subtracting the forward component of velocity of the aircraft from the total 
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velocity vector before driving the cab. The effect can be visualized by the assumption that 
the aircraft flies in a wind tunnel in which the air mass moves at a constant velocity equal 
to that of the forward speed of the aircraft. When the thrust of the engines is adjusted so 
that the aircraft is stationary with respect to the wind tunnel, the situation is aerodynamical ly 
identical to the motion of the aircraft in a windless airmass but with constant forward velocity. 

In -consequence of these considerations, the control system for the motion simulator in 
a situation in which one-to-one motion simulation is possible has the form shown in Figure 
2-3. The aircraft acceleration and angular velocity vectors, in body axes are transformed 
to accelerations of the cab along its axes of travel and to cab gimbal rates, respectively. 

The former are integrated once to provide the velocity components of the cab. After sub- 
tracting the constant forward speed v^ of the aircraft, the velocity components are again 
integrated to produce the position commands for the linear drive system. Concurrently, the 
cab gimbal rates are integrated once to yield the cab gimbal angle commands which are used 
to drive gimbal servos. 

When the task or set of tasks to be simulated is such that one-to-one motion is not per- 
missible, however, the linear and angular drive signals cannot be generated as shown in 
Figure 2-3. Instead, it is necessary to "wash-out" some of the aircraft motion before generating 
the cab servo drive signals. The placement of the wash-out system (or "wash-out circuits") 
is shown in Figure 2-4. . 

In the design of the wash-out circuits, three techniques are conventionally used: scaling, 
high-pass filtering, and "residual tilts". Scaling consists of multiplying each component of the 
vector acceleration or angular velocity by a constant scale factor less than unity. This causes 
an attenuated sensation of motion, but the sensed directions of the vectors in the simulator are 
the same as in the aircraft. High-pass filtering is employed to eliminate the d-c and low- 
frequency components of acceleration which lead to large excursions. To compensate for the 
loss of sustained (i.e. !ow-frequency)accelerations due to high-pass filtering, tesidual tilts 
are sometimes used. The idea here is to use the components of the gravity vector in the 
forward and lateral directions which result when the cab is tilted to simulate the sensation of 
sustained acceleration in these directions. 
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Simulator Control Without Washout 
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Simulator Control Showing Placement of Washout Circuit 








The amount of scaling, filtering and tilting which is employed in a particular sim- 
ulation is currently determined empirically, using a combination rules-of-thumb developed 
out of prior experience and adjustments on the simulation of an aircraft in an actual mission. 

It has been generally recognized that this design technique may not use the capabilities 
of the simulator to the fullest extent: that a more systematic approach might provide greater 
fidelity of motion within the confines of the same volume, or might be useful for a larger 
variety of tasks. The purpose of this investigation is to study the possibility of using the 
techniques of optimum control theory as a method of systematizing the design of improved 
wash-out systems. 

If a single control system is to be used for a broad spectrum of tasks, it would appear 
that the desired performance is essentially nonlinear (with respect to the functional de- 
pendence between the wash-out and the aircraft motion). For those tasks or phases of a 
task in which one-to-one simulation is possible, one-to-one motion should be used. When 
one-to-one simulation is not possible, the minimum amount of wash-out which keeps the 
cab within its permissible motion limits should be used. The basic question in designing 
the control system is: What aspects of the aircraft motion must be followed accurately by 

the simulator, and what aspects can be sacrificed without degrading the realism of the 
simulation? It is generally accepted that one of the most important factors governing the 
realism of a simulation is the kinesthetic sensation of motion. Although there is a continuing 
discussion of what "sensed" quantities really are, a consensus of the opinion is that the linear 
acceleration and the angular velocity, as measured with respect to axes fixed in the pilot 
are the most pertinent factors sensed by human kinesthetic sensory organs. Thus, simply 
stated, the control problem is to endeavor to find a control law which, while keeping the 
cab excursions within specified physical boundaries, minimizes the errors in the motion 
sensations. If a suitable cost or penalty function of the errors in the motion sensations can be 
determined in addition to an analytical model for representing the errors in the motion sen- 
sation, then it is possible to apply the optimum control technique. Since this approach is 
considered the most likely candidate to yield the "best" simulator control law, it was adopted 
for this investigation. 
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During the initial phase of the investigation of applying the quasi-optimum control 
technique to the wash-out circuit design (1968-1969), a relatively simple case of one-degree- 
of-freedom longitudinal motion was considered to establish the feasibility of using the technique. 

In the present phase of the investigation, the general six degree -of-freedom motion 
simulation is considered. In order to treat the six degree-of-freedom problem realistically, 
it was necessary to modify the analytical approach used inthe first phase because the earlier 
approach led to hopelessly complicated calculations in the more realistic problem. The problem 
formulation in the present case is substantially more realistic in that a trade-off between the 
angular and the linear motions is allowed to utilize the effect of "residual tilt" and that the 
washout control system is independent of the aircraft dynamics, although "good" setting of the 
parameters in the wash-out circuit may depend on aircraft dynamics as well as the specific 
mission. 

In the following subsections, the different phases of the analysis leading to the simulator 
control design are discussed in detail and these include the following: 

. Definition of the various coordinate systems and their mutual transformations. 

. Formulation of the problem in a manner suitable for application of the optimum 
control technique including a discussion of the notation and pertinent quantities. 

. Development of various cost or penalty functions of the errors in the sensed 
motion which provide a realistic measure of the "goodness" of the performance 
and, on the other hand, are mathematically tractable as a performance indices 
for the optimum control technique. 

. Application of the quasi-optimum control technique to obtain the simulator 
control law. 

. Discussion of an implimentation scheme for realizing the simulator control 
design. 
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2.2 Coordinate Systems and Transformations 


There are two sets of axes of significance in the motion simulation: a set of "inertial " 
axes fixed with respect to the ground and a set of axes fixed in the vehicle and moving 
: with it: 

(1) . - Inertial Coordinates - The assumed "inertial" reference frame is a cartesian co- 
, ordinate system with its origin at the center of the moving-base simulator track assembly 

and oriented so that each axis coincides with one direction of translational motion of the 
simulator. The positive directions of the x, y, z axes are chosen, respectively, to coincide 
with forward, right side and downward motions. The translational motions of both the cab 
and the aircraft are defined with respect to the same inertial coordinate system. 

(2) Body Coordinates - The body axis forms a cartesian coordinate system fixed with 
respect to the vehicle where the origin is located at the pilot's seat in the cockpit of the 
vehicle. The directions of the axes are in the same sense as the inertial reference (i.e. 
forward, right, and down), but with respect to the pilot rather than ground. Since there are 
two vehicles, the simulator cab and the actual aircraft, there are correspondingly two sets 
of body coordinates. The two sets of body coordinates are illustrated in Figure 2~5 where 
subscripts "A" and " c " are used to denote "aircraft" and "cab", respectively. 

The displacements of the vehicles from the inertial frame are denoted by position vectors 
r^ and r^_ * . The angular orientations of the vehicles are described by the Euler angles 
(yaw, pitch, and roll) relating each respective system of body coordinates to the inertial co- 
ordinates. Any vector defined in the inertial coordinates can always be transformed into 
either of the systems of body coordinates by means of an orthogonal direction cosine matrix. 
Elements in the direction cosine matrix are functions of the Euler angles where the particular 
functions depend on the sequence of rotations from the inertial directions to the directions of 
the body axis that have been adopted. 


* In this report, variables with an arrow (->) on top denote physical vectors with magnitude 
and direction (velocity, force, and etc.) whereas variables with a bar ( - ) on top merely 
denote column vectors with three elements where the three elements may not be the x, y, z 
components of a physical vector. 
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In accordance with the physical structure of the existing gimbal assembly, as 
shown in Figure 2-6, the angular orientation of the cab is described by the three gimbal 
angles : 

Outer gimbal : Pitch ( 0 ) 

Middle gimbal: Yaw(jJ)) 

Inner gimbal: Roll (cp) 

Thus, for computational convenience, the Euler angle transformations in this report for 
both the cab and the aircraft are defined according to the following sequence of rotations: 
(Pitch — Yaw — Roll). In mathematical notation, let 




0 . 

t 


- Euler angle vector 


n^ = a vector in the inertial coordinate system 

= the transformed vector of n. in the body coordinate system 

1/ V 

where £ = A for the aircraft 
t = c for the cab 

then the vector u is obtained from n by the transformation 


fx . ~ C(\. ) n . 

' 7 . 7 7 . 


( 2 . 1 ) 


where the transformation matrix C(X. ) is the orthogonal direction cosine matrix given by 


C( V - 


COS )l) cos 0 

' 7 7 . 


sin i|>. 


sin 0. sin 0. - sin ib. cos m cos 0 cos <t>. cos tb 

v i v i v v 

cos <t> sin 0. + sin O . sin ib . cos 0. - sin <t> . cos ii). 

•f 7 7 . * 7 7 7 *7 


- cos A sin 0. 

y i v 

sin cos 0 + cos 0 sin sin 

tv V V 

cos 0 cos 0 - sin 0. sin tb. sin 

v i v T v 


( 2 . 2 ) 
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Figure 2-6 

Gimbal System for the Ames Simulator 
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(2.3) 


The orthogonal property of C(\ ), i.e. 

tf 

f 

C _1 (x.) -C'(L) 

u u 

is used frequently in subsequent calculations 

It is noted that this particular sequence of rotations is specific to the physical structure 

of the existing gimbal system and is not the generally adopted sequence of rotations used in 

the description of aircraft motion. The conventional sequence of rotations is (yaw, pitch, 

/ 

roll). If the conventional set of Euler angles are denoted by X A = ((p A » 9^/ 4> A )/ then 
these are related to the Euler angles X A = (<P A / 0 A , j p^) , used here as defined in (2.1) 
and (2.2) by 


‘P. - tan 
A 


j { sin <P^ cos - cos <P^ sin 0^ sin ^ 
cos <P^ cos ^ + sin <P£ sin 0£ sin ip* 


0 A - tan 


-1 


sin 0^ 


cos 0* cos ip* 


-1 


Ip A - sin (cos 0* sin jjj*) 


(2.4) 


Let CO be the angular velocity vector of the vehicle with p , q . and r. as its com- 

i v t v 

ponents in forward, right side and downward axis of the body coordinates, respectively . The 
components p. , q. and r. , customarily refered to as "roll rate", "pitch rate" and "yaw rate", 

respectively, are related to the gimbal angle rates X . = (<p . , 0 . , Ip . ) 7 by 

V l jy ■ if 




The superscripts (-1) and (f) used in conjunction with a vector or a matrix denote, 
respectively, the inverse and the transpose of the particular vector or matrix. 
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i. 


sin ib. 

~z> 

COS <P COS Ip. 

V 

■sin <p . cos ip . 

l, ' V 


sin <p. 

1 / 


COS <p. 
7 / 


<p* 


(2.5) 


For riotational convenience, let 

1 

F(X.) = 0 

7 / 


sin if) 

cos <P cos ip 

v z> 


0 -sin <p. cos jl). 

t v 

then (2.5) can be written in vector notation as 




sin (p. 


cos <p. 

. 7/ 


( 2 . 6 ) 


(2.7) 


The matrix F( \ ) is not orthogonal in general since the three components <P , 0 and 

i • 1/1/ 

.I) of the Euler angle rate vector \ . are not orthogonal components of a physical vector. The 

v 

inverse of F(X ) which will be encountered frequently in subsequent calculations, is given by 


F \> * 


1 - tan jJj^cos <p tan j^sin <p^ 

0 cos V* /cos ip^ -sin cp^os ip^ 


si n <P. 

7 / 


' cos c p, 
V 


It should be noted that the orientation of all the angular quantities used in this report are 
defined in accordance with the right-hand convention as shown in Figure 2~6. 
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2.3 Mathematical Formulation of the Control Problem 


The general control problem is formulated as a set of first order differential 
equations relating the state of the dynamic system to the control variables and also an 
explicit expression of the performance index in terms of the state variables and control 
variables defined in the differential equations. The objective of the system design is 
to find the control law for computing the control variables from the state variables 
which minimizes the performance index. The system of differential equations defining 
the dynamic system under investigation is described below. Several -choices for the 
explicit expression of the performance index is given in the next section. 

The first step in the mathematical formulation of the dynamics is to derive the 
differential equations describing the motion of the simulator cab plus explicit expressions 
relating the sensed motion felt by the pilot in the simulator cab to the actual motion of 
the cab. There is a corresponding set of variables and differential equations describing 
the actual and sensed motion of the reference aircraft. The equations for each vehicle 
are distinguished by the subscript A or c which will be used throughout the report to 
indicate reference to the aircraft or to the cab, respectively. 

Since we are interested in how well the cab motion duplicates the aircraft motion, 
it is convenient to introduce a third set of variables defined as the difference between two 
corresponding quantities for the cab and for the aircraft. These variables will be referred ' 
to as the "error" quantities and will be denoted by the absence of a subscript. 

Since we are generally dealing with variables or quantities composed of three com- 
ponents, the development will utilize three -component vectors. In the definitions given 
below, it should be noted that, for both the cab and the aircraft, some of the quantities 
are defined in their respective body coordinates and the other quantities are defined in the 
inertial coordinates fixed with respect to the ground. To express these quantities in a for- 
mat convenient for an optimum design approach, we let 
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= position vector of vehicle in inertial reference axes. 


ix 


ty 


v iz 


=7. = veclocity vector of vehicle in inertial reference axis. 




0 


i 


= Euler angle vector relating vehicle body axes to inertial 
reference axes. 


U £0 



U t 0 

= 


u ui> 



f. 

Z/X 



f. 

ty 

• • 

r i 

- g 


- \ — Euler angle rate 

v 




Oi 3 

1 .x 


a 


ty 


a 


I/Z 




"sensed" specific force 

reading from accelerometers mounted at the pilot's seat 
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*1 = 




0 


iy 


v z 


- a . + g - "unbiased " sensed specific force. 


a). - 


P i 


- F(X.)u - "sensed" angular velocity vector of vehicle. 

2/ u 


- reading of rate gyros mounted on the vehicle's body axes. 


where 


i - A for aircraft 
i = c for cab 


0 

0 

9 


- gravitational acceleration vector in inertial reference 


9 - 32.2 ft/sec 2 .= 9.81 m/sec 2 

C{\ ) = orthogonal transformation matrix which transforms a vector in 

' inertial reference axes to a vector in body axes as defined in (2.2) 


F(X^) = transformation matrix which transforms an Euler angle rate vector 
into an angular velocity vector as defined in (2.6) 


The error quantities which indicate the difference between the motion (both actual 
and sensed) of the cab and the aircraft are defined by 



X = K ' X A 

• • • 

j= V“a = x c" x a =x 
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In the above equations, it is noted that the components of the vectors r , v a. , g 

u 1/ ^ 

are defined with respect to the inertial reference axes whereas the components of a - 0 

V 

are defined with respect to the corresponding body axes. 


The quantity fj . , obtained by adding g to & , is considered a more salient 

is is 

measure of translational acceleration sensed by the human pilot rather than a , since 

l 

it is argued that a pilot who is accustomed to being in a 1 g environment does not really 
sense the effect of a 1 g force acting vertically downward with respect to himself. In 
other words, the pilot's "vertical accelerometer" only senses deviations in the specific 
force from the normal 1 g force. In the subsequent development, the unbiased sensed 

specific force j8, will be considered to be the sensed translational acceleration felt by 

t/ 

the pilot. 


Obviously, the sensed errors j3 and w are a consequence of the actual motion 

• • • 

errors r and X . The actual motion errors are in turn generated by the presence of the wash 
opt signals. In other words, from (2.8), the relation 



X = u 

establishes the fact that 

a = translational acceleration washout signal 
u = angular rate washout signal 


From (2.8) and (2.9), it follows that the differential equation governing the actual 
motion of the cab are 
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( 2 . 10 )/ 


The differential equations (2.9) represent the dynamic system required for the 
optimum design procedure where the washout signals u and a are the control variables. 

The differentia! equations (2.10), which form the basis for the realization of the designed 
washout control system, suggest a general control system configuration as shown in Figure 
2-7. As indicated in Figure 2-7, the command signals generated as a result of the pilot's 
manipulation of the flight controls in the simulator cab are processed by the digital computer 
to compute the corresponding motion of the reference aircraft. The variables defining the 
motion of the aircraft are in turn processed by the washout circuit to compute a and u. 

The washout signals a and u are then combined with the aircraft motion according to 
(2.10) to produce the command signals controlling the motion of the cab. The two 
integrators shown in Figure 2-7 are necessary since the servomechanisms of the cab drive 
system are designed to follow the translational position commands r and the gimbal angle 
commands X c * A more detailed block diagram showing the physical implementation of 
the control system is given In Section 2.6.1 after the equations describing the washout 
circuit have been derived, (See Figure 2-19) 

As discussed previously, a simulation is regarded as perfect if the sensed motion 
errors are zero (j3= to = 0). Consequently, the goal of a wash-out control system is to 
minimize 8 and to • It is noted, however, that since 8 and <*> are physical vectors with 
magnitude and direction, the question of whether a particular value of 8 ( or to ) is 
smaller than another value of 8 (or to ) in the minimization process requires further 
interpretation. This question is resolved by introducing two scaler functions, 

and M (to ), which in some sense, measures the "size" of 8 ond to with respect to 

<0 

the fidelity of the simulation. If these penalty functions are to be used in the performance 
index required for the optimum control technique, then the functions M and M must 
possess the following properties: 


<*> 
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SIMULATOR 



Figure 2-7 

General Configuration of Simulotor Control Syitern 









0M 


M >0, — - i f 0, >0 for l / 0 

P 9 0 ajr 

3 M 0 2 M 

M =0, S = 0, — e >0 for B = 0 

* 3 0 3 f P 


3 2 M 


( 2 . 11 ) 


3 M 

M > 0, *!L 0, 


3 2 M 


oj 


u> 


3 w 

3 M 


M = 0, 


to 


-a - o. 


3 to 


3 2 M 

to 


3 oj 


a -2 
3 oj 


> 0 for oo f 0 


> 0 for to = 0 


(2.12) 


From (2.11) and (2.12), it follows that the penalty functions M and M increase 

monotonically with increases in the sensed motion errors and to , respectively. A more 

detailed discussion of the penalty functions, including explicit expressions for M and 

B 

M is given in the next section, 
to 


Another aspect necessary to complete the formulation of the control problem is 
the consideration of physical constraints. The cab motion is, in general, constrained 
by limits on the translational distances and by limits on the gimbal rotations, both the 
angles and perhaps the angular rates. Nevertheless, experiences with the existing 
simulators has indicated that except for the excursion limits, the other limitations 
seldom cause difficulties. For this reason, and for mathematical convenience, only 
excursion constraints will be considered. 

In order to account in the performance index for the constraints on the cab motion, 

a penalty function L(r ) is defined in a similar fashion to M and M except that 
c p iO 

L(r c ) depends only on the magnitude of r^. The optimum control technique requires 
that the penalty function be chosen to possess the following properties: 
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(2.13 ) 


L ( 'c )s 0 '^T >0 

3T 2 

C 


Several specific forms of L(r c ) will be considered in a subsequent section. 

Based on the penalty functions, introduced above, a performance index S is defined as 
t + T 

$ I (M+kM +6 L) dr - 
t j3 co 


(2. 14 ) 


where t - present time, T - fixed time span, and k, 6 are constant adjustable weighting 
factors. The control laws a and u are determined by means of the optimum control 
technique so that the scaler S is a minimum (or near minimum) for any given values of 
k and £ • For the trivial case, £ = 0, i.e., when no boundary constraint is imposed on 
the cab motion, the control law should provide a perfect simulation. For large £, the cab 
excursions being highly penalized, the control law should generate abnormal motion cues. 
In this sense the parameter £ can be regarded as, a control on. the amount of wash-out, ad- 
justed in accordance with the available distances permitted for travel of the simulator. 


The parameter k which weighs M in the performance index S in (2.14) provides 

co 

a means of adjusting the trade-off between translational and angular motions, since the 

effect of "residual tilt" is embodied in the linear combination of the penalty functions 

M and M . The technique of residual tilt can be illustrated by an example where the 

poo _ 

reference aircraft is flying at a trim attitude (x^ = 0) with a forward acceleration of x^ 

( y'^ = = 0) as illustrated in Figure 2-8a . Suppose that the simulator cab has no trans- 

lational motion, but is pitched downward by an angle 0^ as illustrated in Figure 2-8b. In 
the pitched down condition, the simulator pilot will sense a forward translational acceleration 
of g sin 0^ as a result of gravity. If the pitch angle 0^ is chosen so that x ^ = 9 sin 0 


then the simulator pilot would sense the same motion as felt in the actual aircraft,if it is 
assumed that the simulator pilot is unaware of or ignores the fact that the cab is pitched down 
ward. There is also an upward translational acceleration of g (1 - cos 0 ) sensed by the 
simulator pilot as a result of pitching downward, but this upward translational, acceleration 
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is much smaller than the forward translational acceleration as long as the pitch angle does 
not become too large and thus can be ignored. Therefore, by pitching the cab upward or 
downward, it is possible to simulate a translational acceleration along the longitudinal axis 
(x) of the aircraft. Similarly, by rolling the cab, it is possible to simulate a translational 
acceleration along the lateral axis (y) of the aircraft. Clearly, residual tilt can only sim- 
ulate translational accelerations which are smaller than 1 g (g = 32.2 ft/sec 2 = 9.81 m/sec 2 ). 

The disadvantage of residual tilt is the unrealistic angular motion required to im- 
pliment the technique. In cases where high fidelity of translational motion sensing is de- 
manded, but the limited cab motion imposes a severe handicap, the technique of residual 
tilt may result in a better overall simulation of the sensed motion even though abnormal angular 
motion is imposed in order to compensate for the limited translational accelerations. The 
problem of determining the optimum trade-off between translational and angular cab motion 
which yields the greatest fidelity of the simulation forms the primary basis for using the 
optimum control approach. 




• • 

*A 


g sin Q 

C 

— ► 


— ♦ 


0 " 

0 

0 = 

0 

A 

0 

C 

-9 (1 - cos 0 ) 

L c . 


(a) Reference Aircraft (b) Simulator Cab 


Figure 2-8 

Illustration of the Effect of Residual Tilt 
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To summarize, the equations describing the optimum control problem for de- 
signing the washout control are 


Dynamic System: 
• • 

— * — * 
r = a 

X = u 


Performance Index: 
t + T 

S = f (M +kM +6L)dT 
J $ w 


(2.15) 


(2.16) 


Find : 

a and u so that S is minimized 


Simulator Motion : 



(2. 17 ) 


In solving the optimum control problem, it will be assumed that a flight simulation 

always starts from a trim condition and the cab attitude is aligned so that the initial 
• • 

conditions r(t)=r^(t), r(t) = r^(t), and \(t) = 0 hold. It is also assumed that all 
aircraft quantities needed for the control law can be obtained from the computer sim- 
ulating its motion. 
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2.4 Performance Indices 


One of the major task in attempting to apply an optimum control method to the 
control of the six degree -of-freedom simulator is the search for an appropriate perfor- 
mance index. To start with, the state-of-the-art dictates that a precise definition of 
what constitutes a sensed motion cue is still an open physiological question. Even with 
the assumption that the motion cues sensed by the pilot are the same quantities as can 
be measured by three linear accelerometers and by three rate gyros, the fact that these 
quantities are vector valued makes it difficult, if not impossible, to find a single functional 
expression which provides a realistic measure of how well the simulator is able to duplicate 
the motion sensations. The search for a realistic performance index in functional form is 
further compounded by the necessary requirement of mathematical tractability. As often 
the case in the application of the optimum control technique, a final form of the per- 
formance index is chosen among other promising candidates because of its functional 
simplicity. 

As indicated previously, the basic form assumed for the performance index is 
t +T 

S = f (M + k M +€L) dr 
0 co 

In this section, analytical expressions for the penalty functions M , M and L will be 

p 00 

defined in accordance with (2.11) - (2.1.3), and their properties analyzed. It is emphasized 
that the penalty functions discussed below were selected from a larger set of candidates, 
some of which may appear physically more realistic but were excluded from further con- 
sideration because of their functional complexity. 


(2. 18 ) 
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2 . 4.1 Sensed Acceleration Penalty M 


A 


To establish a meaningful penalty function in terms of the measurable sensed quantitiesi 
3 A and 3^ (Figure 2-9 ), it is observed that the following properties should be considered: 

(a) M should be scaled relative tojS A / i.e. If all vectored quantities 
in M are multiplied by the same constant, the value of M should 

P _ . ,S 

remain unchanged. In other words, fora given error 3 , the penalty 
is higher for small 3 A and lower for large 

(b) M should depend only on the orientation of 3 relative 

p c 

to 3^ anc l not upon the absolute orientation per se. 

(c) It should be possible to represent M as a function of the magnitude of the 

_ P _ _ 

error vector | 3 j and the phase angle £ between 3 and 3^ . 


Based on these propositions, a general form of M which thus would be suitable as 

P 


a penalty measure is 


M =f + K f (f) 

8 m \ 0 A / 0 P * 


( 2 . 19 ) 


where the functionals f and f account, respectively, for the magnitude and phase errors, 
m p 

and the weighting factor K is used to adjust the relative importance between the magnitude 

P 

and phase errors. The effect of the parameter on the penalty imposed by for differences 

between the sensed motion of the cab and the aircraft, can be shown by a vector diagram 

comparing the sensed acceleration vectors 3 and 3 * • An example of such a vector diagram is 

given in Figure 2-10 where the vectors 3 c and 3^ are normalized with respect to the magnitude 3^ of 

— 3^y(3^-=| 3^-|— yand-where the unit vector-3^/3^-is used"as a-basis of-reference~ln"Figifre 

2 - 10 , six different possible values of the normalized vector denoted by the indices I 

through 6 are shown for comparison with 3^/3^. The circles with their center at the tip of the vector 

3 a /3 a depict the contours for which the magnitudes of the error vector 3/3^ remain constant. 

Obviously, if K were zero, the penalty imposed on sensed cab motions represented by the vectors 
P 
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are the same despite the fact that cases 1 and 4 have a phase error of £ , while cases 

5 and 6 have zero phase error. With nonzero K , however, the penalty corresponding to 

p 

1 and 4 is greater than the penalty corresponding to cases 5 and 6 by an amount of K f (£). 

P P 


In Figure 2-10, the cab motions depicted by case 1 and case 4 have the same 
magnitude error and phase error; the same is also true for the pair of cases 2 and 3 
and the pair of cases 5 and 6 . For cases 1 , 2 , 5 , the magnitude of the cab motion 
is smaller than that of the aircraft, whereas the corresponding cases 4, 3, 6 have larger mag- 
nitude of motion as shown from the lengths of the vectors. In practice, however, the 
magnitudesof the cab motion is in general smaller than that of the aircraft. An excessively 
large magnitude of cab motion at a particular time interval will have to be compensated 
by an excessively small or even negative motion thereafter in order that the position of the 
cab does not cross the boundary limits. Thus, it would appear logical to place higher penalties 
on cases 3, 4, 6 where the magnitudes of the sensed cab motion are larger than for the aircraft. 


The above discussion can be summarized by sketching contours of constant penalty, 

M = constant, as shown in Figure 2- 11 . 

P 

Another important property that should be considered in selecting M is that of 

P 

assigning greater penalty for large phase errors. This is related to the notion, that main- 

taining in the same general direction as is at least as significant, if rot more so, 

than magnitude errors. Thus, in the region where the phase error is greater than, say,-^-, 

the contours of M should display a pattern as shown in Figure 2-12. 

P 

With these requirements in mind, various choices of M are discussed below. 

P 
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(A). Penalty Function "A" 


M =-i- 

0 * 2 
Pa 


A 


f- + Kg[i90 A ( 1 -cos|')] 


n = l, 2, 




2 »\ "M/>, 


where 


0 = P £ / /3 = | £ | , etc . 

To aid in plotting the constant M contours, let 

A 

X = cos £ 7 

^A 


= -J_ 

^A 


y = — - sin I' 


then (2.20 ) is rewritten as 


M £ t (x2 + y 2 ) + Kj K 7 + y 2 -x)' 


Figure 2- 13 shows the M contours for n = I, K = 0.02, 0.5 and 12.5 

P 0 


M. 


0 


- -y ( x 2 + y 2 ) + (\Z X 2 + y 2 . x ) 


Figure 2-* 14 shows the case for n = 2 and K =0.5 

0 

^j-Jx 2 ±y ? )_jLK 0 


( 2 . 20 ) 


( 2.21 ) 


( 2 . 22 ) 


(2. 23 ) 


■(2—24—) 
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(B) Penalty Fund-ion "B“ 


r 2 




i= _A_ 

2 4 


+ K ^ M- 


— * _ 

* V 


Let 


then 



cos £ ' = 




sin £ ' 




]_ 

2 


(x 2 + 


y 2 ) 


+ K, 




(2.25 ) 


(2. 26 ) 


(2.27 ) 


and the set of curves for n = 1 and =0.5 is shown jn Figure 2-15 


(C) Penalty Function "C" 


M = -£ 


ft' 


2 ^ 


2 + K e 


.fc ‘ K > 4 “«') 



T ( * 2 + y 2 > + V 2 a " K i x)! 


(2.28 ) 


(2.29 ) 


37 



Figure 2- 16 shows the M contours for different combinations of K and K . 

p pi' 


(D) Penalty Function "D" 


- 'V' -V 

* 2/Sa 


'* K ,( A '-*I 


= _ji 


,+K 4’- K i -V 

- 


(2. 30 ) 


In terms of x and y. 


% = t (x 2 + y 2)n +K ^ 0 “ K i x)2] 


(2.31 ) 


Figure 2- 17 shows the M contours for different combinations of K and K. . 

p pi 


2.4*2 Sensed Angular Penalty M 


to 


It is assumed that the angular penalty M takes the same form as the acc'eleration 

penalty M . Therefore, all of the discussions concerning M also applies to M . The 
P < _ P ' W 

various possible choices used for the functi.on are summarized below: 


M •= + K 

to 0 2 a) \ to 

2 W A 


W A w 


to, 


n = l,2. 


(2.32 ) 


to 


M = + K |1- 


: to. 


to 


^A 

W c W A 


n 1,2, ... 


(2.33 ) 
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Figure 2-16 

Contours of Constant M for Penalty 

P 

Function "C" 


Figure 2-17 

Contours of Constant for Penalty 
Function "D" 










(2.34 ) 



(2.35 ) 


2.4.3 Excursion Penalty U r c ) 


The basic requirement in the choice of L(r c ) is that of assigning heavy penalty 
for large cab excursions and small penalty when travelling well within the allowable 
boundary. There are two possible types of the limiting penalty described below. 

(a) "Soft" Limiting Penalty; 

2n 

+ I 


1 


L(r c )= i 




, n a 1 


(2.36 ) 


where d^, d^, d^ are adjustable parameters, whose values depend on the actual lengths 
of the allowable excursion limits. For ri = 1, (2. 36 ) reduces to 


with 


Kr>= * r. d7 c 


D = 


x 

0 


0 

1 


(2.37) 


0 


0 



(b) "Hard" Limiting Penalty: 


where 


L(7 ) = L +L + L 
c x y z 



|x c l<d x 

K 1 2d x 



Figure 2- 18 illustrates these two types of penalty functions. 


(2.38) 


(2.39) 
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j 2.5 QuasirOptimum Washout System Design 

| In this section, the optimum control problem, as formulated in Sections 2.3 and 

2.4, will be solved using the quasi-optimum control technique. Although it would be de- 
I sirable to consider all the penalty functions defined in Section 2.4 so that the relative merits 

j of the resulting washout systems could be evaluated by actual flight simulation, only one 

case will be considered in detail due to the limited amount of time available. In particular, 
the penalty functions considered in this report are (2.20), (2.32) and (2.36) with n = 1, i.e. 



In subsection 2.5.1 , the maximum principle is applied to reduce the optimum control 
problem to a two-point boundary value problem whose physical interpretation is discussed in 
subsection 2.5.2. Finally, in subsection 2.5.3, the quasi-optimum control technique is 
applied to solve the two-point boundary value problem for the controls a and u in closed form. 


2.5.1 The Two-Point Boundary Value Problem 

The~control _ probienrto be^treat'e'd“c6risists _ 6f"the error dynamics' 


r - a 


X - u 


(2.43 ) 
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and the performance index 


S = ^ (M + k M + £ L) d T (2. 44) 

P W 

which is rewritten in state variable form, by letting r Q = S 

r n = M + k M + £ L 
0/3 u> 

• 

— ♦ — ♦ 
r = v 

v = a (2. 45 ) 

• 

A - u 

r = 1 

€ = 0 

The initial time is taken as the present time t and the time interval T is assumed fixed. 

The assumed initial conditions are 

r o ( 0 =o 

r(t) =r A (o 

v(t) = v A (t) (2. 46 ) 

X(t)= 0 
T (t) = t 
€(t) = € 

and the assumed terminal conditions are 

r 0 (t+T)=f ree , X(t + T) = free 

7( t + T) = free , r (t + T) = t + T (2. 47 ) 

v( t + T) = free , € (t + T) = free 


44 



The Hamiltonian for the dynamic system described by (2.45 ) is 


h = p n (M + k M + € L ) + p' v + p' a + j5' u + p 
0 £ cu r r r \ t 


(2.48 ) 


where the adjoint variables p are defined by 


v° 


3L(? ) 


p r * " €P 0 ~W 


P v - 


*"Pr 


3M (o,X,t) 3M (u,X, T ) 


A 


+ k 


CO 


3 X 


3 X 


(2.49 ) 


P r - 


3M (a, X, t) 

Hr +k 


3 M (u, X,r) 3 L ( r, T ) 


CO 


3 T 


+ 6 . 


3T 


p 6 


subject to the boundary conditions 

P 0 (t + T) = -1 

P r (t+T) = 0 

P v (t+T) = 0 

P x (t+T) = 0 

p (t + T)= free 
T 

P £ (t +T) = 0 


(2.50 ) 


The solution of the two-point boundary value problem defined by (2. 45 ) — (2. 50 ) 
is highly complicated in its original form due to the irrational form of the penalty functions 
(2. 40 ) and (2.41). It will be shown that by means of a nonlinear transformation, the original 
problem can be reduced to one which may be more readily solved and which in turn provides 
physical insight into the structure of the system. 
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The first step is to apply the maximum principle to obtain the optimum controls 
a and u. This entails computing the partial derivatives of the hamiltonian- h, given 
by (2. 48 ) with respect to a and u , respectively, and then setting the derivatives to 
zero. The resulting equations are 


3 M 


-P. 


da 

3 M 


co _ - 


3 u 


= ? X A 


(2. 51 ) 


(2. 52 ) 


where the adjoint variable p^ was easily found to be p^ = - 1 from (2.49) and (2.50). 

Substitution of M and M from (2. 40 ) and (2. 41 ) into (2. 5? ) and (2. 52 ), re- 

jo co 

spectively, results in an expression which can be solved for c? and u. The final equations 
for a* and u (see Appendix I for the details of the derivation) are 
K 


a = 


1 - 


JL 




«i p v + 2 K e %C ( X c )p v *K l 




(2. 53 ) 


u - 


1 - 


CO 




2 c-J/r 


■1 - ,2 


p' F"'(X )[F'(\ )]"' p,A" +2K 3/ [ F'(X )]"’ p. A + K 


A ^X 


co 


CO 


-1 


co 2 A F _1 (x )[F'(X )] ’p A +K F _1 (X ) co. 

A c c X co c A 


+ f ’ I( V“a' j a 


(2.54 ) 


or, in functional form 

a -a(p v , X , A a ( T )' °A^ T ^ 
u = u(p x / X , X a (t), 0 a (T)) 
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(2.55 ) 



Using the optimum control laws (2. 53 ) and (2. 54 ), the closed -loop 
two-point boundary value problem can be obtained by substituting these controls into ' 

(2. 45 ) and (2. 49 ). It is noted from (2. 45 ), (2. 49 ) and (2. 55 ) that the pertinent 
variables are T, v*, \, p^, p^ and p^ and that they are uncoupled from the rest of the state 
and adjoint variables. In order to obtain the closed-loop equations in terms of the pertinent 
variables, we evaluate the partial derivatives in (2. 49 ) by noting that 


9M 


3 X 
9 M 


M. 


CO _ 


3 X 


3 X 


9 M 





3 P 

9 M 


co 


3w 


(2. 56 ) 


But, from (2. 51 ) and (2. 52 ), we have 



9 M 


to _ 


3 u 


9 M 


Bl = p 

V 
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a -\' 3M 

liiL) — = p A 

35 1 9co * 


(2. 57 ) 


or 


9 M 


3 0. 

9M 




co _ 


3 oo 



(2. 58 ) 


and- 


= C(X ) 

35* C 


(2.59 ) 


^ = F(X ) 

35 
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ax 

9 to 

ax 


ax 


a f ( x ) 

c - 
u 


(2. 60 ) 


ax 


Substitution of (2. 58 ) — (2.60 ) into (2. 56 ) yields 
9 M 


ax 

9 M 


8 _ 

- a 

c 


3 X 


co _ - , 

u 

c 


ac(x ) 

c 

ax 


a f (x ) 

c 


ax 


C(x c ) p v 


[F'(x c )]~ ] p x A 


(2.61 ) 


Further substitution of (2. 61 ) into (2. 49 ) obtains the closed-loop two-point boundary 

value problem for the pertinent variables 


r = v 
T =a 
X = u 


p - € 

r 


(2.62 ) 


3 L(7) 

9 T 


r r 


(2.63 ) 


P X = °c 


a f (x ) 

c 


ax 




3 c( X c ) 

ax 


. c(x Jp 


where o’ and u are given by (2. 53 ) and (2. 54 ), respectively, and where the initial 
conditions and boundary conditions are given, respectively, in (2. 46 ) and (2. 50 ). The 

objective is to solve for p and p in terms of the state vectors r^ v* and X from the coupled 

v r X 

differential equations (2. 62 ) and (2. 63 ). Substitution of the expressions for p and p 
into the control laws for a and u shown in (2. 53 ) and (2. 54 ), will result in the feedback 
configuration necessary to realize the washout circuit. 
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The nonlinear transformation required to solve the two-point boundary value 

problem involves the introduction of a new adjoint vector y and a new angular washout 

vector n that are related to p and u by 

A 


-i 

y=[F'(\ c )l p x 


n - F(\ ) u 


(2. 64 ) 
(2. 65 ) 


Multiplying both sides of (2.64) by F'(x ) and differentiating with respect to time, we have 


P X = F'(X e )y + u' 


dF'(\) 


3 X 


c - 

V 


( 2.66 ) 


Equate (2.66) to the last equation in (2.63), gives 


-1 dC'(\) 


y = [F'(X )1 a' C(X )p 

c c 0 r c v 


(2.67 ) 


After performing the sequence of matrix multiplications in (2. 67 ), the differential 
equation for the new adjoint vector simplifies to 

y=-a c x C(X c )p v 

= ^(X c )[° c x P v l (2 . 68 ) 


where the symbol "x" denotes the "cross product" of two vectors. For mathematical con- 
venience, (2. 68 ) can be rewritten as (see Appendix I) 


y = ~a c x n 

(2.69 ) 

where 



— 1 1 

f ^ “A + % ?A 


0A | 

* 


(2.70 ) 
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The two-point boundary value problem given by (2. 62 ) and (2. 63 ), is next 
rewritten in terms of the new variables to give 


X = F V )n 
c 


n t) = r A (t) 

v (t)= 7 A (t) 

X (t) = 0 


(2.71 ) 


P r -6 


8L(f ) 

c 


p = - p 
v r 


P r (t+T) = 0 


p v (t + T) - 0 


(2.72 ) 


V = ~ C{X c )l °c X K ] 


y(t +T) = 0 


where the translational control a* is shown in (2. 53 ) and the transformed angular control 


n is rewritten from (2. 54 ) as 


n = 1 


V4y 2 A 2 +2K w ^A + K* 

+ jr A - f(x c ) u A 


W A + K co “A 


(2.73 ) 


2.5.2 Physical Implication of the Mathematical Formulation 

As a result of'the theoretical developments in the preceeding subsection, the optimal 
control problem was reduced to a problem of solving a set of differential equations given 
by (2. 71 ) and (2. 72 ) for the adjoint vectors p^, p^ and y in the terms of the state 
vectors r", v and X • Before proceeding further with the theoretical solution of the problem, 
we will pause at this point to examine the physical .implications of the equations (2. 71 ) 

and (2. 72 ). 
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First, for reasons of simplicity, we will consider the'"linear" case in which 

no "phase" penalty is imposed, i.e. K = K =0. For this case, it follows from 

p 00 

(2 . 53 ), (2 . 70 ) and (2. 73 ), that the system equations given by (2. 71 ) and 
(2. 72 ) reduce to 


* - V/ C 'IV?A- J A 

w 2 

* = ~r r]{l Jv +F " 1 ( ^c )jr A _a A 

3 L(?*) 



“ 4- (a c x a A ) 


(2.74 ) 


(2.75 ) 


Also from the linearity assumption, the closed-loop cab dynamics (2.17) reduces to 



jAs noted' previously, the parameter k is designed to account for the relative im- . 
portance between the angularand the translational motion errors. Thus, it is of interest 
to see how the cab motion would behave for extreme values of k. 
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( i ) k -> m : 


This case heavily weights errors in angular motion which for the limiting situation 
gives that the angular motion of the cab will be identical to the aircraft. This is indeed 
the case, since for k = ®° , (2. 77 ) reduces to 

V f X>=a 

which, upon substitution of the assumed .initial condition X c (t) = X^(t), further reduces to 


X = 
c 


X> 


u>. 



The translational motion equation (2. 76 


), with X = X . / can be written as 
c A 


= «a K+°a + ? 

2 _ *4 

= *A P v + r A 


(2.78 ) 


(2.79 ) 


which shows that the cab acceleration differs from the aircraft acceleration by the amount 

2 w 2 

0 A P v • T he wash-out signal 0 A P y which can be obtained from (2. 75 ) is dependent on 
the excursion penalty L(r^) and its adjustable weighting € . 

( ti )k ->0: 

Direct visualization of the effect of k = 0 is not obvious from the closed-loop 
equations since the terms involving k tend to co as k-*0. Referring to the performance 
index, given by (2.44), it is intuitively obvious that, since no penalty is imposed on the 
angular motion error for k = 0, the cab should assume whatever angular motion is necessary 
to minimize the translational motion error. It is derived in the AppendixlH that this is indeed 
the case. In fact, as k-0, the optimal trajectory tends to approach a "singular" subarc for 
which the adjoint variable y = 0 and for which 
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(2. 80 


«c X *A ~ C{X c ) °c Xot A 


= 0 


The implication of (2.80) is that along the singular trajectory the cab attitude 

X should be so maintained to keep the sensed translational accelerations & and n. 
c . c A 

co-linear for any specific force a^ being exerted on the cab. The fact that the two 
vectors are co-linear provides that the cab motion is either completely in phase, or com- 
pletely out of phase relative to the aircraft. The possibility of the motion being completely 

out of phase arises because K was set to zero and no phase error penalty was imposed. 

P 

Substitution of \ , obtained from the solution of (2.80), into (2.76) yields the 
c 

corresponding translational motion for the cab. 

Thus, for intermediate values of k, the cab motion will always have some translational 
and angular motion errors, and it is a subject of experimental study to determine what 
particular values of k will give an acceptable compromise between the two extremes. 


2.5.3 Quasi-Optimum Solution of the Two-Point Boundary Value Problem 
Simplified Control: 

In applying the quasi-optimum control technique, the system is first approximated 
by a simpler model for which the exact optimum solution can be expressed in closed form. 
This simplified control law is then corrected to account for the difference between the 
original system and its simplified model. A convenient simplified system in the present 
i cation js obtai ned for € = O.^i .e . -byJgnoring-the-physical constraint— For 
the two-point boundary value problem, given by (2. 71 ) and (2. 72 ) reduces to 



53 



= V 


= a 


£, -F-’IXJS. 


r (t) = ^ 

v 

s ^ 

x 5 (t ) -0 
5 (t+T) = 0 


= 0 


rs 


rs 


= 0 


Pv S ' Pr ‘ 


p <> + n 

r vs 

; (t+T) = o 


~ = -a xC(X J P vs 

CS . the sW^ed system. 

, V h as been added to s.gn.fy 

where the subscrip 


Solution 


of (2. 82 ) immediately yields 


= 0 


rs 


= 0 


vs 


r s '° 


73 


co ^' y .'-" 0 - 53 ,OTdl2 ' 

? =C!k c5 )? A - a A 

S -S A - FU ^“ A 


. . /p Pi ) and noti ng 
c (0 84 ) ,nto ^ ’ ° 

Substitution of U - 84 1 


), th. Aimplified controls 


are 


(2.81 ) 


(2.82 ) 


(2. 83 ) 


(2.84 ) 


«x.=x s ^ A .-^’ hec '“ ed ' 


CS 


loop error state equations 


(2. 85 


b =C'(X e + X A )«A '°A 


t =F-'a. + A A )ffA' c A 
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It is noteworthy to observe that the simplified controls in (2.84) and the correspond- 
ing error states in (2.85) are not zero, even though no excursion limits are imposed on the 
cab motion. However, the "sensed" errors are zero which can be seen by substituting 
(2.84) into the definitions of 0 and go given by (2.8) . 

P = *cs~«A = C(I cs ) K + °V ~ “ A = 0 

3 = *cs ' S A = F(I cs )(u "s +G A ) ' "a = 0 

■ V. . 

An examination of the error state equations (2.85) reveals that the nonzero simplified 
controls and the nonzero error states are a consequence of a nonzero initial attitude error 
X $ * If it happens that is initially zero, then the solution of the last equation in (2.85) 
yields \ $ (t ) = 0 which upon substitution into the remaining equations in (2.85) and in 
(2.84) results in the errors and the controls being zero. 

In summary, when no boundary limitations are considered (simplified control problem), 
a one-to-one simulation, i.e. 0 = go = 0, can be achieved for arbitrary initial conditions 
by using the wash-out controls in (2. 84 )• 

The solution of the simplified problem is required to calculate the quasi-optimum 
control law to be derived in the next section. Unfortunately, there are two major dif- 
ficulties that arise when attempting to obtain analytical solutions to the simplified problem: 

(A) The transformation matrices C and F , as defined in (2. 2 ) 

and (2-6 ) are nonlinear. 

(B) In order to solve the two-point boundary value problem, the 
quantities n a and X^ describing the aircraft motion 
must not only be known for the current time t but for time in 
the future of t. 
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To overcome these difficulties it will be assumed that the attitude error X between 

s 

the cab and the aircraft is sufficiently small to permit the approximation 


C 'W’ f A S 

F ’ 1(X cs )5 A F ~ 1( *A )jS A ‘ u "a = ° 


( 2 . 86 ) 


Such an approximation is not unreasonable since, with no excursion limits imposed on the 
simplified system the attitude error will tend to decrease from its initial value. If the initial 
attitude error is zero, then it will remain at zero. 


From (2. 86 ), (2.85) can easily be solved to give 
% (T) =r s ( t )+v s (t) (r - t) 

v* s (t) = ?(t) (2.87 ) 

X, (r) =X s (t) 


Quasi-Optimum Control: 

The quasi -optimum controls are obtained by first deriving the correction factors 
l|) , ip and jp so that the adjoint variables p^ r p^ and y of the original system can 
be approximated by 

p = p +1 = i 
r rs T r r r 

p = p + i = i (2. 

r v vs T v r v 

y = y + ip = ip 

Next, p and y are substituted into the control laws for a* and n given by (2. 53 ) and 
r v 

(2. 73 ), respectively. 
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I CL 


and y are given by 


It is shown in Appendix H that the quasi-optimum solutions for the costate variables 


P = - € 
v 


t + T t+T T 

T J ^r^-J J— dXdr 


t dr 
t+T 

y S y(0-J (a c XT])d T 


t t 3 r 


(2.89 ) 
(2.90 ) 


In order to evaluate the integral in (2.89), where the analytical expression for the 
function L( r) is given in (2.42), knowledge of aircraft's position r^(r) for time T in the 
future of the current time t(r £ t) is required. To resolve this problem, a weighted Taylor's 
series expansion of r^(r) is used to extrapolate the future position of the aircraft from its 
current value r^(t). The resulting expansion is 


where 


r A (r) = r* A (t)MT-.)r>) + ifL -p;w 

max 


(2.91 ) 


&A max ' s maximum value of 0 ^ that must be determined (approximately) a priori. 
Utilizing (2. 91 ), gives 


9L 


9 r 


€ =0 


= D 


€ =0 


= o ( r + r A ) 


€ - 0 


“ D(f ' s +f A ) 


D r* 
s 


(t) + ( T - t) r* (t) + rt (t) + (t - t) rT (t) + 


(t - tf 


A' 1 ' 2 

2 |ST 

__A_max- 


*>>] 


d r 

c 


+ k - 


cs (.) ♦ (r - f ) r s („ 2 

^ 8 A 

A max 




(2. 92 ) 


57 



Therefore, dropping the subscript "s" we have from (2. 89 ) 



The evaluation of (2. 90 ) will be done numerically during the simulation because of 
the complicated expression for the integrand . 
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2.6 Summary and System Implementation 


In this section the simulator control system designed previously using the quasi- 
optimum control technique is summarized and its implementation is discussed. A block 
diagram showing how the required computations can be implemented during an actual 
flight simulation is given. A FORTRAN IV subroutine which performs the computations 
shown in the block diagram is also described. 

2.6.1 Summary of the Simulator Control System Design 

Recalling from Figure 2-2, the complete six degree -of-freedom simulator consists 
of the simulator cab, aircraft computer simulation, simulator control system, and cab 
drive system. The flight commands resulting from the pilot's manipulation of the cab con- 
trols are fed into the computer simulation of the aircraft. The computed aircraft motion, 
in particular, the translational acceleration r^ and the angular velocity V are the 
inputs to the simulator control system which in turn generates the drive command signals 
.used to control the motion of the cab. In the problem formulation it was assumed that 
the cab drive system consists of perfect servos so that the outputs from the simulator control 
system are the cab position r^_ and the cab attitude \ . The basic input-output equations 

for the simulator control system are 

(2.94) 

(2.95) 


» • » » 



"translation " 
"rotation" 


where 


y t 


z J 

v 


x',= ffmfV 

t t V r 7/ 


i = c,A 

as illustrated in Figure 2-4. Thus, the simulator control system has a total of 6 inputs and 
6 outputs. 


59 


The washout signals a and u are included in (2.94) and (2.95), respectively, in order 
to compensate for the limits on the cab excursions. The purpose of this study is to design the 
control laws for computing a and u . The performance index used in section 2.5 for 
determining a and u was 


-1 


t+T 


(M +kM + £ L)d T 
8 oi 


(2.96) 


where = p 2 / 2^ + ft //9^) 


P 


- & /2a)^ + ^ W A ” w// W A ^ 


L = 

2 c c 


D = Dlag{l/d*,1/d* 1/d*} 

x y z 

This performance index minimizes the sensed motion errors while penalizing large 

cab excursions. This performance index gives rise to 8 scaler parameters T, k, f, K , K , 

8 co 

d^, d^, d^. A large portion of the computer simulation study described subsequently is 
concerned with an investigation of the effects of using different values of these adjustable 
parameters in order to arrive at a range of values that may optimize the performance of the 
moving base simulator. 

In Section 2.5, the quasi-optimum technique was used to design a washout system 
that minimizes the performance index in (2.96). The resulting equations for computing 
a and u are summarized below. 


Translational Washout Signals : 

a = 0 - K g 8 A /0l + <=' (I )5 a - a A 

C = 4F v+ k s c'(x c )S a 


(2.97) 
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Rotational Washout Signals: 

u = F 1 (x c )n 

n = (1 - K^ a j A /o)a+ W A - f(x c )x a 

O = ( oi i A / K ) y + K ^ A 

Adjoint Variables 

F v = - e D((T 2 /2)7 c + dVsjV + 


(2.98) 


y ~ ~a c * rj 


(2.99) 


v =n/d-y A / 5 ) ? |]; A + (y/ A )^ 


Sensed Motion Variables 
« • 

o A = r - g 


«A" C(X A )a A 


8 a -« a + 9 


Transformation Matrices 
f cos ii cos 0 


a c = r c -g 


o - C ^ r ) a c 

c c c 


(2.100) 


a> A F(x a ) x A 


= [ 0 , 0 , g ] (g = 32.2 ft/sec 2 = 9.81 m/sec 2 ) 


cos i}) cos 0 sin 0 - cos 0 sin 0 

C (x) = sin 0 sin 9 - sin 0 cos 0 cos 0 cos 0 cos 0 sin 0 cos 0 + cos 0 sin 0sin 0 
cos 0 sin 0 + sin 0 sin ijj cos 0 - sin 0 cos 0 cos 0 cos 0 - sin 0 sin 0sin 0 

■■ — ; (2rl 01 ) 



f 

1 

sin 0 

\ 

0 

, F _, (I) = 

f 

1 

- tan 0 cos 0 

tan 0s?n 0 

= 

0 

COS 0 COS (J) 

sin 0 

0 

cos 0 sec 0 

- sin 0 sec ^ 


0 

- sin 0 cos il) 

cos 0 

> 


0 

sin 0 

cos 0 


-/ = [ 0 , 6 , 0 ] 
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A block diagram of the simulator control system showing the implementation of (2.94)/ 
(2.9 5), and (2.97) - (2.101) is given in Figure 2-19. Each arrow or path in the block diagram 
corresponds to a 3-component vector quantity. A total of 3x5=15 integrators are required 
to realize the simulator control system. The boxes in the block diagram represent gains 
multiplying the various vector quantities with the one exception of a single vector cross product 
indicated by a bold face X . The gains are composed of both scaler and matrix multiplication 
of the vector inputs to the boxes. The boldface letters C, D, and F indicate matrix operations. 
The matrix D is a constant, diagonal matrix (see (2.96)). The matrices C and F denote co- 
ordinate transformations which are a function of either \ ^ or \ ; the particular choice is 

indicated in the block diagram. Thus, the matrices C and F are nonlinear gains which must 
be computed continuously. To simplify the block diagram, arrows connecting the integrator 
outputs \ and \ to the corresponding boxes containing C and F are not shown. 

Of the scaler gains, some are linear, constant gains while the others are nonlinear, 
time-varying gains which must be computed at each iteration. The nonlinear gains require 
computation of ft ^ which are the magnitudes of the vectors g^, ^ , £/ O, 

respectively. These magnitudes are obtained by computing at each iteration the square root 
of the vector dot product. Again for simplicity purposes, the computation of these magnitudes 
is not explicitly shown in the block diagram. 

In summary, the block diagram shows that the washout signals a and u are a function 
of 

• • • • 

a = a ^A'^A' V r c'VV v) 

' ••• .• 

u = U ^A' ^A' *A' V f ' *c ,y ^ 

where 

X^,7 7^, \ c , y : state variables, of the simulator control system (outputs of the 

integrators). 

• • •_ 

r A' ^A : ' n P uts *° the simulator control system. 

r , X : outputs from the simulator control system, 
c c 
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Figure 2-19 

Block Diagram of the Simulator Control System 
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2.6.2 Analysis of the Linear Design Case 


The general characteristics of the simulator control system shown in Figure 2-19 are 
of interest and would be valuable in selecting the optimum values of-the parameters. Unfortu- 
nately, a detailed analytical analysis is very difficult to perform because of the nonlinear 
cross-coupling between the six degrees- of-freedom. Consequently, the selection of the 
parameters must be performed by computer simulation. However, for the special case of the 
quadratic performance index to be discussed below, it is possible analytically to gain some 
insight into the characteristics of the simulator control system. 

IF K = K = 0 in the performance index shown in (2.96) then the equations for the 
8 a) 

simulator control system given by (2.94), (2.95), (2.97) - (2.101) can.be simplified. In 
this case it can be shown that the equations for computing r and X from 7 ^ and X A are 


V A .V A r^-VA +[C'rX c )C(x A ) - j]a A 


t c = F ''S c ) 


2 /»t— * -* 

F-r- W A / a c X *A , 1 

f £a )x a-it J ~ T ~ dT \ 

L Jo p A J 


( 2 . 102 ) 

(2.103) 


where matrix coefficients in the first equation describing the translational motion of the cab 
are given by 


'1 




1 2 2 

A 2=J eT S A° 

A 3'H« T4 ’A 0/ 4™x 


The second equation (2. 103) for the angular motion of the cab can be rewritten as 
a> A j —J o a T (2.104) 


cu - a r 


2 f- - 
“A| a c X “A, 

= — r~ dr 

0 (*A 
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If £ = 0 then = A ^ 13 ^3 = ® anc ^ ^ rom (2.100) it can be shown that (2.102) 
reduces to oj. - S A or J . 0. Substituting S c • S A in ,o (2.104) immediately gives 

that a) c = 6 ^ or a) 2 0. In other words, if no penalty is imposed on the cab excursions 
then the simulator control system will cause the cab to move so that there is no error in the 
sensed motion. For £>0 there will be errors in the sensed motion. 

v 

An examination of (2.102) shows that the translation motion of the cab is governed by a 

• • • • 

second order system with two forcing terms r^ and a^ g. The second forcing term 

on the right hand side of ( 2 . 102 ) makes adjustments in the translational cab motion due to 
differences in the cab and aircraft attitudes. If \ = then the second forcing term 

C M 

vanishes. The second forcing term occurs because the objective is to minimize the sensed 
errors rather than the actual motion errors. 

• « 

Ignoring the second forcing term in (2.102), the computation of 7_ from is almost a second 

order linear system with constant coefficients; it would be a timei-invariant, linear system except 
2 

for the term in the matrix coefficients Assuming for the sake of argument 

that fj ^ is a constant, then the transfer function from 7^ to 7 is a high pass filter. To 
illustrate, suppose the characteristic equation of the left hand side of ( 2 . 102 ) has real roots 

and then the Bode plot of the frequency response would appear as shown in Figure2-20. 

Figure 2-20 
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The high pass filtering characteristics that evolved from minimizing the performance index 
are intuitively reasonable since one of the common approaches is to attenuate the low frequency 
portionof the aircraft motion while passing the high frequencies. Another feature incorporated 
in the translational equation (2.102) is the presence of the scaling factor (r“A ) multiplying 
the aircraft acceleration. Equation (2. 102) shows that the magnitude of the cab position 
will tend to vary linearly with difference between jand A^. Thus, the adjustment of the 
constant parameter e, T, D and (3 ^ max' that affect the values of A^ , A ^ and A^, in turn 
amounts to the trade-off between the high pass characteristics and the degree of scaling. 

The rotational motion of cab governed by (2.103) or (2.104) encompasses the idea of 
residual tilt. An examination of (2.104) shows that as long as the sensed translational motion 
vectors for the cab and aircraft remain colinear (i.e., x = 0) then the cab attitude 

identically follows the aircraft attitude. As the phase error between rx c and ot ^ deviates 

from zero, (2. 104) induces an angular error in order to help reduce the error ^ and a in 
the sensed translational motion. The trade off between angular and translational errors in the 
sensed motion is governed by the parameter k . As k increases, the level of angular 
error decreases according to (2. 104). 

2.6.3 Computer Simulation Program 

In order to study the performance of the simulator control system design, a FORTRAN 
computer program was developed to simulate the complete six degree-of-freedom motion 

simulator. The simulated time histories shown in Section 3 were obtained by use of this 
program. The computations required to implement the simulator control system design shown 
in Figure 2-19 were performed in a subroutine entitled WASHFL. 

The inputs to the WASHFL subroutine consists of the aircraft motion given by , to , 

and r ^ , the integration step size £t , and an index variable governing the mode of 

operation. The outputs from WASHFL are the cab motion command signals given by 7 and 
™ c 

\ c . The adjustable parameters in the performance index may be entered as input data to 

the program or stored internally in the subroutine in DATA statements. The latter storage option 

permits WASHFL to be used on the NASA AMES computer facility in order to generate the drive 

commands during an actual flight simulation. 
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A more detailed description of the WASfiFL subroutine, including a program listing and 
description of all computer variables, is given in Appendix IV. Two other subroutines required 
by WASHFL to perform the numerical integration of the state variables are also given in '• . 
Appendix IV. 


3. PERFORMANCE SIMULATION 

3.1 Description of the Reference Aircraft Motions 

To have a realistic assessment of the performance of the washout circuit designed in 
the preceding sections it ultimately would be necessary to conduct a full-fledged pilot simu- 
lationusing the simulation facilities at the Ames Research Center. As a preliminary step 
toward such a simulation, digital computer simulation studies were conducted to evaluate 
possible performance of the quasi-optimum washout circuit and to determine suitable ranges 
of the values for the adjustable parameters in the washout circuit. 

The first phase of the computer simulation study was to generate examples of the aircraft 
motion that result in several distinctive "tasks') which the simulator cab attempts to duplicate, 
during typical operations of the motion simulator. The second phase, discussed in the next 
section, was to examine the cab motion resulting from the application of the washout circuit 
for the reference aircraft motion generated in the first phase. 

The dynamics used to compute the reference aircraft motion are for a medium sized 
twin jet transport whose aerodynamic characteristics are summarized in Table 3.1. 

To generate the reference aircraft motion, a general six degree -of- freedom aircraft 
motion simulation program (SIXDOF) was used. This program was developed at the Kearfott 
Research Center and has been used extensively for a variety of studies. The SIXDOF program 
incorporates an autopilot, also developed by the Kearfott Research Center as part of an 
automatic landing study. This autopilot was used to simulate the behavior of a human 
pilot. To use the SIXDOF program, a nominal trajectory for a given typical flight task is 
first defined in the manner required by the program. The autopilot then computes the aircraft 
control -surface deflections required to minimize the difference between the actual aircraft 
flight path and that of the defined nominal trajectory at each instant after the inception of 
simulation run. From the control-surface deflections, the reference aircraft motion is generated 
by the model of the aircraft dynamics. 

The three flight tasks considered in this study are the following: 

Task 1 Tracking Maneuver - This task was created for the purpose of simulating such 
flight operations as mid-air refueling, formation flying, etc. It is also a useful task forevalu- 
ation of quantitative pilot followup error. The reference trajectory of the aircraft and the 
time histories of 24 dynamic variables of the aircraft are shown in Fig. 3-1 for 30 seconds 
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1 

6 e / 6 q / 6. = elevator, aileron and rudder deflections 

It is seen from these trajectories that the flight is quite hectic. This is intentional 
and designed to give the washout circuit a good workout. 

Task 2 Approach Landing Maneuver - Aircraft landing operation is another typical ex- 
ample of simulator application. To generate a reference aircraft motion, it is assumed that 
during a landing approach, at about 50 seconds before touchdown, the pilot "suddenly" realizes 
that the aircraft is 200 ft* (61 m) and the right of a prespecified nominal landing tra- 
jectory and that, instead of aborting the landing as he would normally do under these circum- 
stances, he tries to complete the landing. In view of the relatively large course error, just 
50 seconds prior to touchdown, a rather severe maneuver is called for. The autopilot that 
simulates the pilot performance completes the desired landing successfully 55 seconds after the 
maneuver is initiated. ' 

Figure 3-2 shows the resulting reference aircraft trajectory for those dynamic variables 
listed in the preceding section. 

Task 3 Emergency Pull-up Maneuver - The last simulation mission selected was an 
emergency pull-up maneuver during a final landing approach. In this mission the aircraft was 
assumed to be initially flying the ILS beam. At an altitude of about 50 ft (15.2m) above the runway 
the pilot decided to abort the landing and excuted the pull-up maneuver, in which the aircraft 
pitched up sharply while simultaneously turning away to the right side of the runway. 

A particular feature of this flight operation is the large yaw and pitch angles it induces, 
and, as a consequence, the large deviation in forward (x) acceleration. Figure 3-3 shows the 
resulting trajectories. 


of flight time. These variables are: 

x, y, z = inertial x, y, z excurisions 

x, y, z = inertial x, y, z velocities 

= inertial x, y, z accelerations 

x, y, z 

' (v , a / a = sensed x, y, z specific forces 
x y z 

"■ to, ft -0 = Euler angles 

Ip, Q,‘tp = Euler angle rates 
p, q, r = angular rates 
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TABLE 3.1 

CHARACTERISTICS OF SIMULATED AIRCRAFT 
(a) Physical Properties 


c 

11.08 ft 

(3.377 m) 

m 

125000.slug - ft 2 (169477. kg-m 2 ) 

m 

71.2 

ft 

(21.70 m) 

EH 

120312.slug - ft 2 (163121. kg-m 2 ) 

s 

690. 

ft 

(210. m) 

■ 

234375.slug - ft 2 (317770. kg-m 2 ) 

m 

777.5 

slug (1 1347. kg) 

i 

xz 

8125.slug - ft 2 ( 11016. kg-m 2 ) 

■ 

236.4 

ft/sec 

(72.05m/sec) 

HI 1 

30 ft (9.14m) ahead of 
vehicle CG 


(b) Stability Derivatives 
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Figure 3-3 

Reference Aircraft Trajectory - Emergency Pull-up Maneuver 
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3 .2 Simulated Time Histories of Cab Motion 

The washout control law as summarized in (2.97) - (2.101) contains nine adjustable 
parameters/ namely: 


d 

x 



£ - adjust the amount of linear translational magnitude washout 
k - adjust the amount of linear angular magnitude washout 
- adjust the amount of nonlinear translational phase washout 

K - adjust the amount of nonlinear angular phase washout 
00 

T - adjust the "filtering" characteristics between high-cut-off-and-low-scaling 
and low-cut-off-and- high-scaling 

d - adjust the weighting of £ in x, y, z directions 


a. - limits the maximum scaling 

PA max 

The task of choosing a suitable combination of these parameters is highly complicated in view 
of the multidimensional performance requirements and the interrelationship between these 
parameters and performance. For the purpose of facilitating a systematic determination of 
the performance and the adjustable parameters, the following sequence of cases were con- 
sidered during the course of simulation study. 


Case 1 Linear translational washout only (K - K =0, k-»co) 

|3 w 

Case 2 Linear trc.islational and angular washout ( K = K =0) 

jS O! 

Case 3 Nonlinear translational and linear angular washout (K =0) 

u> 

Case 4 Nonlinear translational and angular washout 

It should be noted that these cases were studied in consecutive order and a parameter which was 

found suitable for a previous case was not altered for the following case unless it was necessary 

to maintain the cab excursions to stay within confinement. In other words, the simulation is 

not exhaustive for every case due to limitation in available computer time. It is conceivable 

that better results than those shown in this section may be achieved by a more exhaustive 

simulation. 


Of the three tasks considered, the tracking and landing maneuvers received a fairly 
thorough investigation but only Case 1 and Case 2 were considered for the emergency pull-up 
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maneuver. Typical time histories of the cab motion and that of corresponding aircraft motion 
for these tasks and cases are shown in the following pages and their corresponding values of 
the parameters and figure numbers are tabulated in Tables 3. 2-3. 4. 


In the Figures 3-4-3-29, the cab motion is represented by solid lines and the corresponding 
aircraft motion, by solid lines with circles. The dynamic quantities shown in these figures are 


x,y, z = inertial x, 
p x , p z = "unbiased "x 


y, z excursions 
, y, z sensed specific forces 


<p, ft, ih ~ gimbal angles 
p, q, r = angular rates 

The cab excursions shown in these figures are all confined to within ± 10.0 ft. (3.05m) and 
aircraft excursions, which are in the order of 1000ft (300m), were purposely "trimmed" down 
from the actual values shown in Figs. 3-1-3-3 in order to permit plotting both cab and air- 
craft excursions by the same scale. Let x r . , , z,. be the aircraft excursions shown 

fig 'fig fig 

in Fig. 3 _ 4-3-29 and x . y , z , be the actual aircraft excuisions as shown in 
act 'act act 

Fig. 3 -1-3-3, then for 


Tracking task: 


Landing task: 


Pull-up task: 


*cct 

x fig 

+ 30.0 +236. 7 1 (ft) 

^ X fig 

+ 9.14 + 72.15 

Km)] 

^act 

= y n< 

3 




Z act 

Zf '9 

- 5001.5 (ft) 

[ Zp 

f'g 

- 1524.5 (m)] 


*act 

a> 

u- 

X 

1! 

+ 30.0 +243. Ot (ft) 

[ x r . 
f'g 

+ 9.14 + 74.0t 

(m) ] 

^act 

= y fig 

+ 200.0 (ft) 

[y ng 

+ 61.96 

(m) ] 

Z act 

= z fig 

-718.8 (ft) 

[z fig 

-219.1 

(m) ] 

X act 

= x fig 

+7320.5 + 120.0 1 (ft) 

Cx r r 

f'g 

+ 2231 .3-+ 36.. 

58 t (m) ] - 

^act 

= y fig 

+ 90.0t (ft) 

[y fig 

+ 27.43 (m) 1 


Z act 

z fig 

- 127.2 - 13. 5t (ft) 

Cz ng 

- 38,77 - 4.1 1 i 

t (m)] 
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Values of Adjustable Parameters for the Tracking Task 





































Figure 3-4 

Relative Motions of the Aircraft and the Cob - Tracking Maneuver 
- Case 1 (linear translational washout only), - T = 500,0 
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Figure 3-5 

Relative Motions of the Aircraft and the Cab - Tracking Maneuver 
- Case 1 (linear translational washout only). - T = 40.0 
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Figure 3-9 

Relative Motions of the Aircraft and the Cab - Tracking Maneuver 
- Case 2 (linear translational and angular Washout). - T = 2.0 







Figure 3-10 

Relative Motions of the Aircraft and the Cab - Tracking Maneuver 
- Case 3 (nonlinear translational and linear angular washout). - T= 500.0 
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Figure 3-14 

Relative Motions of the Aircraft and the Cab - Tracking Maneuver 
- Case 4 (nonlinear translational and angular washout). T = 40.0 
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Figure 3-18 

Relative Motions of the Aircraft and the Cab - Landing Maneuver 
- Case 1 (linear translatio'nal washout only). T = 2.0 
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Figure 3-1? 

Relative Motions of the Aircraft and the Cab - Landing Maneuver 
- Case 2 (linear translational and angular washout). T= 30.0 
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Figure 3-25 

Relative Motions of the Aircraft and the Cab - Landing Maneuver 
- Case 4 (nonlinear translational and angular washout). > T = 30.0 
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Figure 3-26 

Relative Motions of the Aircraft and the Cab. - Landing Maneuver 
- Case 4 (nonlinear translational and angular washout). T= 10.0 
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3.3 Summary and Discussion of Simulation Results 

2 

Of the nine adjustable parameters £, k, K , K , T, d^, d , d and , the three 

p (i) y ^ max 

parameters k, K , K were used to define the four cases in Section 3.2, the parameters 

P M 2.i 

€, d , d , d were used as direct control of the cab excursions, g. , which limits the 
x y' z “A max 

the magnitude of maximum scaling, were kept constant throughout the simulation study 
after initial trial runs, whereas the remaining parameter T was used to alter the nonlinear 
filtering characteristics of washout circuit for the translational motion. 

It is noted from Section 2.6.2 that for large T, the system tends to behave as a high- 
pass filter with low cut-off frequency and high scaling of the input aircraft accelerations, 
while for small T, the system is more like a high-pass filter with higher cut-off frequency 
and low scaling effect. In general, the low-cut-off-and-high-scaling effect is reflected in 
the sensed cab motion with good phase and poor amplitude relations with that of sensed air- 
craft motion. On the other hand, the high-cut-off-and-low-sCaling effect tends to give poor 
phase but with better amplitude relations. These effects are visibly evident for all the tasks 
and cases considered in the present study as will be pointed out in the sequel. 

Another important property for assessing performance is the "onset" characteristic, i.e., 
the ability of the cab to follow the initial aircraft translational acceleration. The onset 
motion generally occurs in the vertical (z) direction since the turns are normally well co- 
ordinated and forward speed variation is more often kept at low level. A glance at the 
part (6) of the Figs 3 -4-3-29 reveals excellent onset following characteristics of the washout 
system for all tasks. 

As for the angular motion, it is generally agreed that the deviation between the sensed 
angular motion (i.e., angular rates) of the cab and the aircraft should be preferably kept small 
since the human motion sensor tends to be more sensitive to the deviation in the angular motion 
than to the deviation in the translational motion. In the present washout system, the amount of 
angular motion deviation is controlled by the adjustment of parameter k, which has been kept 
at a value that keeps angular motion error small. 

Thus, the major effort in the simulation study is to obtain a range of the adjustable 
parameters to provide suitable choice of translational motion with limited motion deviation. 
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In the discussion that follows, attention will be directed to the response of the sensed 
translational motion 0^, and ^ in part (4), (5) and (6)', respectively, of Fig. 3-4- 
3-29. 


Tracking Task : 

Three values of the parameter T were consideredin each case to illustrate the trade- 
off between th e good-phase-poor-amplitude and poor-phase-good-amplitude responses. Take 
Case I (Figs. 3-4, 3-5 and 3-6) for example: at the initiation of the simulation, all cases 
show excellent onset following; after I sec, however, small T (Fig. 3-6 (6)) is’seen to cause the 
sensed cab motion p^ to go in an opposite direction to that of sensed aircraft motion 
(phase error) while for large T (Fig. 3-4 (6)), p^ remains at the same direction as that of 

• On f be other hand, at about 4 seconds, large T causes p to reach a magnitude 
of 8 ft/sec/sec whereas, for small T, the peak p reaches only 4 ft/sec/sec. Again - 
during the following period of 4 sec, large T. provides, better amplitude, response with the 
cost of phase error while small T provides better phase relation with the cost of larger 
amplitude attenuation. An intermediate value of T (Fig. 3-5 (6)) results in a compromise 
between the two extremes. 


The reduced peak amplitude response caused by large T , which may not be favorable in 

z-direction, is definitely favorable in y and x-directions. Because of well coordinated 

> • • . 

turns, the aircraft has a small lateral sensed acceleration , but f be cab motion which can- 
not be coordinated within the available maneuvering region shows some spurious lateral 
acceleration p . If is seen that the peak p for large T (Fig. 3-4 (5)) is significantly 
smaller than that for small T (Fig. 3-6 (5)). The same effect is also true for the forward 
acceleration p as can be seen from Fig. 3-4 '(4'), 3-5 (4), and 3-6 (4). 

The effect of the adjustment of T on the cab motion discussed above for Case j 
(Figs. 3-4, 3-5, 3-6) can be extended to cover Case 2 (Figs. 3-7, 3-8, 3-9), Case 3 
(Figs. 3-10, 3-11, 3-12), and Case 4 (3-13, 3-14, 3-15). 


The introduction of angular washout by the adjustment of k in Case 2 is designed to 
utilize the residual tilt for the purpose of improving translational motion. This effect can 
be visualized by comparing the corresponding time-histories for Case 1 (Figs. 3-4, 3-5, 3-6) 
and Case 2 (Figs. 3-7, 3-8, 3-9). For the limited amount of washout provided, a general 
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improvement in x and y directions are visible from the figures. 

In Case 3 (Figs. 3-10, 3-11, 3-12) an additional nonlinear washout signal was pro- 
vided for the purpose of improving phase relations in the translational motion. The effect is 

relatively moderate for the combination of parameters used in the simulation, however, the 
nonlinear washout does provide an alternative to the linear cases. 

In addition to the translational nonlinear washout in Case 3, a nonlinear angular 
washout signal was utilized in Case 4 (Figs. 3-13, 3-14, 3-15) for the purpose of improving 
angular phase relations. Since the angular phase relations are nearly perfect in all cases, no 
significant improvement can be expected from the nonlinear signal. The figures are presented 
here for the sake of completeness. 

Landing Task : 

The response characteristics of the washout system discussed in the foregoing paragraphs 
for the tracking task also applies to this task in general. A major difference between the 
tracking task and the landing task is that the former has a relatively symmetric vertical 
acceleration whereas the latter is asymmetric as can be seen from part 6 of the figures. For 
a motion with asymmetric acceleration profile, the high-cut-off, low-scaling characteristics 
of the filter can be utilized more advantageously to eliminate the d-c component of the 
acceleration. Thus, substantially smaller values of T can be used in the present task than 
had to be used in the previous task. 


Another consideration in the determination of suitable combinations of the adjustable 

parameters is that because of the larger vertical excursion than the lateral excursion in 

the landing task, the weighting on the amount of translational washout (d , d , d ) need 

x y z 

not be equal in all three directions as was the case in the tracking task. In a direction in 
which large excursion occurs a larger weighting is needed to restrain the cab from moving 
beyond the physical boundary, and in a direction with small aircraft excursion the weighting 
can be smaller to permit the cab to move in a larger volume. 


Again, three values of T were considered in each of the four cases and its effect to 
the cab response is similar to that in the tracking task as can be observed by comparing the 
three figures in each case: Figs. 3-16. 3-17, 3—18 for Case 1, Figs. 3-19, 3-20, 3-21 for 
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Case 2, Figs. 3-22, 3-23, 3-24 for Case 3 and Figs 3-25, 3-26, 3-27 for Case 4. Good 
vertical onset following and improved lateral and forward performance are visible in all 
cases. 

Pull-up Task : 1 

Substantially different aircraft motion, reflected in the presense of large forward accel- 
eration variation, in this task from the preceding tasks provides a good test of the adaptability 
of the washout system. Only Case 1 and Case 2 were considered and only one value of 
T was simulated for both cases. 

Referring to Figs. 3-28 and 3-29, a perfect following of the aircraft motion by the cab 
for the first 6 sec, in which the aircraft was descending at a trim condition, is as expected. 
After the commencement of the pull-up maneuver, the sensed cab accelerations display a good 
onset following before subsequent deviation from that of aircraft motion. The effect of high- 
pass characteristics is most visible after 20 sec: the high frequency component of the aircraft 
vertical acceleration was faithfully preserved by the cab motion, while d-c component of 
the acceleration was effectively removed. The effect of residual tilt can be seen by comparing 
Figs. 3-28 and 3-29. A general improvement in cab acceleration is visible in Fig. 3-29 
particularly during the period beyond 20 sec. 
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4. CONCLUSION AND RECOMMENDATION 


This investigation is concerned with the design of a six degree -of- freedom motion 
simulator control system. The problem considered herein is a generalization of the problem 
considered in an earlier study [ 1] in which control systems for a two degree-of-freedom 
motion simulator were designed by applying the quasi-optimum control technique described 
in [2.3] . Although the same general method was employed in the present investigation, the 
mathematical formulation of the physical problem here is substantially more general and more 
convenient, in that, (a) the reference aircraft dynamic was not necessary in the formulation, 
(b) the effects of residual tilt and phase error were incorporated in the nonlinear performance 
i ndex . 

The results of this investigation may be interpreted as demonstrating that: 

• The quasi-optimum washout control system can be effectively used for a 
wide variety of flight simulation tasks, 

• For each simulation task, the cab excursions can be readily confined to 
within any specified value and the characteristics of motion sensation can 
be varied by adjustment of constant parameters. 


It is believed, however, that the adaptability of the washout system has not been 
fully explored in the simulation study reported here. The general presence of spurious 
lateral sensed acceleration p throughout the simulations, for instance, may 
be an inherent limitation of motion simulation in a finite volume and may have to be 
present irrespective of the manner in which the washout signals are generated. On the 
other hand, the possibility of improvement through better combination of adjustable param- 
eters cannot be discounted. Another option, which is logically promising but has not been 
simulated in this study, is the employment of different filtering charactistics of the washout 
system in each direction. It was pointed out in Section 3.3 that, to achieve an appropriate 
washout effect, the filtering characteristics embodied in the washout system should be altered 
by adjustment of parameters from task to task. But for any one task, the filtering characteristics 
were maintained the same in all x, y, and z directions despite the fact that the acceleration pro- 
files are significantly different in each direction. Thus it is conceivable that tailoring the 
filtering characteristics in each direction in accordance with each acceleration component may 
ultimately provide improved performance. 
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Based on a comparison of the results obtained in the present study and those achieved for 
the two degree-of-freedom case which has received favorable comments from pilots who made 
actual flight tests, it is not unreasonable to conclude that we have achieved a design which 
may be preferable to those achieved by use of conventional washout techniques. 

In view of these results the following further effort can be recommended: 

• More exhaustive simulation study of the washout system design reported here. 

• Experimental assessment of the washout system by means of actual flight 
simulation at the Ames Research Center. 

• Extension of present design approach to include the consideration of 
human kinesthetic sensor models. 

The importance of experimental evaluation with actual pilots cannot be minimized. 

It is possible to scrutinize an unending number of time histories without knowing for certain 
whether one design is better than another, because the characteristics that are being sought 
in the time-history have not been pinpointed. Thus, actual pilot experiments are ultimately 
indispensible for complete evaluation of washout system performance. 
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APPENDIX I 


Derivation of Optimum Controls a, u 

In this Appendix, analytical expressions for the optimum controls a* and u will be 
obtained from (2.51) and (2.52), respectively. 
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2 “ \ “A 


“A “ 


GO > 


ff-4) 


Since the solution of u from (1-2) and (1-4) is essentially the same as the solution 
of a* from d~l) and d~3), we will first solve (1-1) and (1-3) for a*. 

Applying chain rule to (1-1), we have 


3 a 



3M^ 

9^ 


= c'U c ) 






a-s) 


where use of the following relation was made 



Equating the right hand sides of (1-1) and (I -5) gives 



+ I S 




= c (‘Vv 


a-7) 


Next substituting (t-6) into (1—7) and solving for a*, yields 


Pa P 


.2 - 




— 2 +c '(X c )a A -o- 

K /3^A 


a-s) 


To eliminate the (3 term in (T-8), (1-7) is rewritten as 


K 


K 


+ -JL\« = -A 

Pa P ] aP 


P ~ 2 Pa + C{ \ ) K • 


a-9) 


p 


Next, compute the inner products of the vectors on either side of equal sign in (1—9); the re- 
sulting scaler equation may then be solved for £ to obtain 


P ={*, 


+ 

P 12 2 


C( K>K + K e 


a- io) 


Substitution of (1—7 0) into ([-8) gives the desired result 

K. 


a - 


1 


P 


2 P 2 v+ 2K s r A 'ca c )p v *< 2 j 


(^ Pv + K S C ''\>? a ) + C '< X c ) “ A -°1 


a-n) 


Following the same procedure results in a similar expression for the angular wash-out, given below 

K 


u - 


1 - 


-Ml. 


\K ^ ^(^[^(^c)!" 1 + 2K aj ^[F'(X c )]- 1 p x A - K 


F' 1 (X c ) [F'XX.)]" 1 p x /k +K w F- 1 (X c )5i 


+ f_1 (X c Ko a - u A 
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APPENDIX H 


Derivation of Quasi -Optimum Correction Factors 

The optimum control problem defined in (2. 94 


) — (2. 97 ) is summarized below: 


r_ = + k M + 6 L 

0 j3 co 

/ 

o_^ 

-f 

II 

o 

> 

r Q (t + T) = free 


"V 

n 

<i 


/ 

f*(t) = r A (t) , 

• F*(t + T) = free 


— ♦ _ 

v - a 


f 

v>(t) =^ A (t), 

v(\ + T) = free 

(n-i) 

vl- 

ii 

C| 


/ 

A (») = 0 , 

X(t +T) = free 


T = 1 


/ 

r (t) =t , 

T (t + T) = t + T 


€ =0 


/ 

€(t) = 6 , 

€ (t + T) = free 


h=p 0 (M 9 

+ KM +6 L) +p' v 
co r 

+ P V ' 3* +P X ' a 

+ P r 

ai-2) 

o 

II 

o 

■Q- 




; P 0 (t + T) = -1 


P r = -€P 0 

9 L 
dr* 



; p f (t + T) = 0 


P v = - P r 




; P y (t + T) = 0 


[ 

— l'[F' 
9 X J 

7 " i^/ 

1 , P X °C 


' p x (t + ^ = 0 

(H-3) 

P T = -P 0 1 

r 9M S 

2- +k 

< 3T 

3 M 

— ) 
9 T / 

; p ( t + T) = free 
T 

; p, (t + T ) = 0 



Let the vectors x and p represent all the state variables and the corresponding 
adjoint variables, respectively, 

I 
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ar-4) 


The corresponding quantities for the 


"simplified" problem (6 = 0) are denoted by 


the subscript "s": 


' r 0s" 


' P 0s' 



___ 

r 


p 

s 


r rs 


__ 


V 

s 

P s = 

P vs 

X c 



s 


P Xs 

T 



S 


P rs 

_0 


0 


(n-5) 


In this appendix, we will derive a correction matrix M so that the original adjoint 
vector p is given by 


P - P c + ! P 


= p +M £ 


ai-6) 


where jjj is the correction to adjoint vector p and £ is the correction to state vector x. 


x = x + F 
s * 


ax-7) 
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The matrix M is conveniently partitioned as follows: 


M = 


00 

m „ 
Or 

m 0 v 

m ox 

m 0 T 

m oe 

Vo 

M 

rr 

M 

r v 

M 

rX 

m 

rr 

*r€ 

v0 

M 

v r 

M 

V V 

“vX 

m 

vr 


x 0 

m n 

\r 

M 

X v 

M xx 

"Xt 

™X€ 

Vo 

m 

T r 

m 

T v 

"tX 

m 

TT 

"V* 

€0 

m, 

€ r 

m 

€ v 

m €X 

m, 

€ T 

m €€ 


(H-8) 


The elements denoted by capital M are n x n matrices, and m's denote vectors. The elements 
appearing in (H-8) — not all of which are required for quasi-optimum control law — are to 
be found with the aid of the auxiliary equations for i p and £ 


£ = H £ +H jl) 

’ xp PP 


lb — H £ - H ib 

Y xx ■» px r 


(n-9) 


or, by use of the matrix Riccati equation: 


d M/d r = M H + H M + MH M+H 

xp px pp XX 


(It- 10) 


The coefficient matrices H , H , H and H in dI-9) and (11-10) are obtained by 

xp px pp xx 

partial differentiation of the Hamiltonian h in (II— 2 ) with respect to x and p where the 
resulting derivatives are evaluated at x^ and p^ . The first partials of h are: 
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OI-il) 



ar-i2) 



In order to evaluate the second partial matrices about the "simplified" solution, it 
is noted that for the simplified problem we have 


6 = 0 

l =<? = £= 0 

V« + °" a = c ' ( V“a 

V G+G A* F " ,(I c )jr A 

P 0 (r) = -1 

p r (T) =0 
P v (t) = 0 

P x ( r ) = 0 

p (t) = p(t + T) 

T 

p At) ~ J L(r(r))dT 

^ L 


or- 1 3 ) 


ar-14) 


CEL-15) 


r 0 (r) = 0 

r*(T ) = r*(t) + vT(t) (r - t) 
v(r)=7(t) 

X(r) =I(t) 

t(t) = r 
€ (r) ~ € 


(n-i6) 
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where the subscript "s" has been omitted in (Et-13) — flt - 16). In the subsequent 
development, we will be using only the “simplified" solution, and hence for con 
venienc the subscript "s" will be omitted unless otherwise indicated. Thus, 


j-] = — j} = |_j/ _ 

PX " 0x0p € =0 xp ‘ 


0 

0 

0 

. 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

I 

0 

0 

0 

0 

0 

0 

-o' 

c 

f0C 
1 9 ^. 

f 

c 

— / 
-u 
c 

3_F' 
9 1 

[Fl" 

0 

0 

0 

0 


3 a*" 
dT_ 


r 3 u 

L 9T . 

t 

0 

0 

L 

— 

0 

0 

0 

0 

0 


C-17) 


'4 





















[5J3 




0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

do 

a p 

r v 

0 

0 

0 

0 

0 U 

0 

0 

0 

0 

0 

0 

0 

0 

H = a2 , h 
xx dx 2 € = o 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

a 

a 

0 

0 

a 













































































Substitution of QI"17) — (E - 19) into the auxiliary equations 
results in 




^r Cv 



(a) 

(b) 


1°. t ■ + 
8 T € f 

9a - 

a - K 

9p v 

(c) 

+ lit 
9 r V 

+ —h 

9p x 

-(d) 


(e) 

(f) 




«€ 



(a) 

(b) 

(c) 



ar-20) 


ar-2D 
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with boundary conditions 

t o (t+T) = "Po (t + T)d(t+T)=0 

^ r (t +T) = - p r (t +T)d (t +T) =0 
^ v (t + D = - p v (t +T)d(t +T) = 0 
4^(t + T) = - p^(t + T)d(t + T) = 0 

* € (t+T) = -p 6 (t +T)d(t + T) = 0 

i (t +T) = -f(t +T)d(t +T) =-d(t +T) 
’r 

and p ' (t + T) £ (t + T) = x '(t + T) (t + T) or 
l/> (t + T) = 0 

T 


(a) 

(b) 

(c) 

(d) 

(e) 

(f) 


01 - 22 ) 


(g) 


Since the optimum control laws a and u in (2. 53 ) and (2. 54 ) are only a function 

p^ and p^ , respectively, the elements in the M matrix (H-8) of primary interest are 

those affecting the correction of p and 5 . The correction factors il) and il) are 

v_ r \ 

written in terms of the elements of M from (H-6) to give 


if) = m ? ar-23) 

i b = M n £ + M £ £ + M 7 + m t +m t (EL-24) 

*v vO 0 vr r vv s v vX ’X vT s r v€ 

if).=M n t n + M + M £ +M X f +m t +m t (H.-25) 

X0 Xr Xv^v XX ^X Xr X€ *6 v ' 


Those elements of the matrix M , appearing in dI-24) and (U-25) will be evaluated 
as follows: 
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From (H~21a ), (H- 22a ) and 01” 23 ), and the fact that M is symmetric/ it is 
evident that the first row and column of M vanish. Hence: 


*n =G 


(E-26 ) 


m v6 ~ m xo 


From (IT- 21b ) and the fact that = £ , we have: 

</) (r) = j/j (t) + € J ~ d x 

r r t 3r 

and hence, from (H-2lc ), we have: 


(n-27 ) 


r Mol 

4 > v (r) = ^v (f) “^r (f)r " € I I 

t t 3r 

Using boundary conditions (H- 22b ) and (H- 22c ) in QX- 27 ) and (H- 28 ) yields: 


(n-28 ) 


J t+T ~ dX = 

t 3 r 


and 


4) v (t) = 


* +T a! * +T t 31 

T j — dX-J J — d X d T 

t 3 r* t t 3 r* 


(n-29 ) 


= m f 
v£ 


A comparison of (H— 29 ) and (II- 24 ) indicates that. 


M = M =M=M =0 
v r vv v X v T 


Therefore, from (H- 6 ) and (II- 15 ), we have the quasi-optimum solution for p 
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t + T 



J T -^d-Xd ; r 

3 r , t t 3 r* 


CEL- 30 ) 


To obtain ip , it is observed from (H-2id ) that the differential equation govern- 

ing tp is exactly the same as the original equation for p in (IT- 3 ), which, as was 

A A 

shown in Section 2.5.1, can be simplified by the transformation defined in (2- 64 )• 
Therefore, proceed by letting 
-1 _ 

r q = [ F'( x c ) ] l p x (n- 31 ) 


then, from (2. 66 ) - (2. 68 ), (H- 2 1 d ) is transformed into 

♦ . , 

y =-C(X Maxi] (H- 32 ) 

q c c ~ r 

and integration of (II- 32 ) yields 

t 

ai-33 ) 

it immediately follows that 

y ( T ) = y (T) = y (0 - J C(\ ) [a* x p ] d r 01- 34 ) 

» j. c c V 

The quasi-optimum approximations of p^ and y in (H- 30 ) and (H- 34 ) respectively 
are the necessary terms in deriving the quasi-optimum control law in Section 2.5.3. 


y_( T )=y n (0- J C (\ ) [ a* x |ji ] d t 

q q c c T v 


Since from (H- 6 ) and (II- 15 ), we have the original p as 

A 

3 *x 
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APPENDIX m 


Singular Optimum Control 


Consider the optimum control problem described by (2. 45 ) - (2. 50 ). If 
k = 0, then the Hamiltonian becomes 

h =P 0 (M^ + eL) + p' v* + p y a*+p^ c + P T (m- 1 ) 

Using (2. 61 ), the adjoint equation (2. 49 ) can now be rewritten as 


P 0 = -1 


P r -6 


3L 

ar 


p = - p 
v r r 


d M 


P X = 




3X 


3 c( x c ) 

3 X 


c (*c>P v 


9M 


BL+c ± L 


3T 


3r 


P, -L 


(nr- 2 ) 
an- 3 > 
an- 4 ) 

an- 5 ) 

an- 6 ) 
an- z ) 


The optimum controls a* and u are obtained by computing the partial derivatives of the 
Hamiltonian h with respect to a* and u and equating the derivatives to zero. 


3M 

+ p = 0 an- 8 ) 

35 

p x = o an- 9 ) 
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Substitution of an- 9 ) into (til- 5 ) results in the relation 


a' 

c 


ac{\ ) 

c 


3 X 


C(x )p =0 

c r v 


an- io) 


which can be rewritten as 
a' C(\ ) 


3C(X ) 

C 


d\ 


C(x c )p v = o 


or, using (2. 71 ) with K = 0, it becomes 

P 

3 C(x V* 1 

c 


«c C ^o J 


3 X 


a A =0 


an- ii ) 


After performing the matrix multiplication, (DI-.ll ) can be reduced to the component 


form 


a , a A “ at ~ 0 


l cy> 2 u cz u Ay 
^cy “Ax ' a cx a Ay } s! " ^c C0S *c + (a cz“Ax " a cx a Az } COS ^c C “ 
+ ^cy a Az ~ a cz a Ay^ sin ^c = ° . 


an- 12 ) 


U a* ~a a. )cos0 +(a ax ~a a* )sin<t> =0 
’“ex Ay cy Ax c cz J Ax cx Az c 


Since the set of equations in Q3Q- 12 ) should be satisfied for all values of <t> c and ^ , 

it follows that 

a cx a Ay- a & y“Ax = ° 

a cy a Az" a cz“Ay =0 (BL- 13 ) 

“cz a Ax" a cx a Az = ° 
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which implies that 

a c x ot A = 0 (DI- 14 ) 

l 

It is noted that for the singular case, the control a is solved from (EH- 8 ) in 

exactly the same way as for the nonsingular case. The control 0, however, can not 

be obtained by a direct application of the maximum principle. The equation (El- 9 ) 

obtained by the maximum principle results in an indeterminate value of u. To obtain 

u for the singular case, first solve the set of equations (HI— 12) for the cab attitude 

0,0 and ih . Next, differentiate the expressions for 0,0, and A and use the 
c c T c c c ~ c 

definition (2. 8 ) to obtain u . 

It is noted that the cross product relation (HI- 14 ) can also be obtained by- 

i ' 

neglecting the angular error dynamics and treating the angle error X as a control variable . 
In this case, the system dynamics (2. 45 ) become 


(DI- 15 ) 

T - ' 


€ =0 



ception of the motion is obtained without regard to angular perception. 
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APPENDIX IV 


FORTRAN SUBROUTINE OF THE SIMULATOR CONTROL SYSTEM 

A FORTRAN subroutine called WASHFL was developed for implementing the simulator 
. control system design shown in Figure 2-19 and given by (2.94), (2.95), (2.97) - (2.101). 

This subroutine was used in the computer simulation study described in Section 3 and also 
has thej capability of being used on the computer facility at NASA-AMES to generate the actual 
drive commands for the moving base simulator. 

this appendix contains program listings of the WASHFL subroutine and two other sub- 
routines, GINTR and DRPRC, required by WASHFL to perform the numerical integration of 
state variables of the simulator control system. In addition to the program listings, a brief 
description of how to use the WASHFL subroutine and a description of all computer variables 
are given. 

SUBROUTINE WASHFL 

Implement the simulator control system design. 

Procedure; See the block diagram of the simulator control system design in Figure 2-19. 
Usage: CALL WASHFL (IMODE, H, APHI, ATHT, APSI, AXFM, AYFM, AZFM, 

PA, QA, RA, CPX, CPY, CPZ, CPHI, CTHT, CPSI) 

where the inputs are 

index controlling the mode of operation; 

= 0, bypass; = -1 , initialization; = + 1 , simulation. 

integration step size (sec). 

aircraft roll angle (rad). 

aircraft pitch angle (rad). 

aircraft yaw angle (rad). 

aircraft inertial acceleration in the x direction (ft/sec ). 

aircraft inertial acceleration in the y direction (ft/sec z ). 

aircraft inertial acceleration in the z direction (ft/sec ). 

aircraft roll rate (rad/sec), 
aircraft pitch rate (rad/sec), 
aircraft yaw rate (rad/sec). 


IMODE 



APHI 
ATHT 
APSI 
AXFM 
AYFM 
; . / AZFM 

. PA 
QA 
RA 




and where the output's are 


CPX 

cab inertial displacement, command in the 

x direction (ft). 

CPY 

cab inertial displacement command in the 

y direction 

(ft). 

CPZ 

cab inertial displacement command in the 

z direction 

(ft). 

CPHI 

cab roll angle command (rad). 



CTHT 

cab pitch, angle command (rad). 

■V 


CPSI 

cab yaw angle command (rad). 




Remarks 

The input aircraft Euler angles APHI, ATHT, and APSI are defined according to the 
conventional sequence of rotations: yaw-pitch-roll. Since the simulator control system was 
designed using the sequence' of rotations: pitch-yaw-roll, a transformation is included in 
WASHFL to convert the aircraft Euler angles from the former sequence to the latter sequence. 
The FORTRAN statements containing the transformation are designated by the word MAYBE 
in columns 74-78. If the aircraft Euler angles defined according to latter sequence are 
available in the computer simulation of the aircraft then, by deleting the FORTRAN state- 
ments identified by MAYBE in columns 74-78, those Euler angles can be used in the calling 
argument. tl 

The output cab Euler angles are defined according to the rotational sequence: pitch- 
yaw-rol I . 

An improved first-order Euler method is used to perform the numerical integration. An 
integration step size of H = 0.05 sec was used in computer simulation studies. 

Subroutines required : GINTR 

A list of all computer variables used in the WASHFL subroutine is given below where the 
variables are partitioned according to scaler or array type. 

Scaler Computer Variables 

Description 

Double precision value of T. 

Double precision value of H. 

t + T is the terminal time in the performance index. 


Computer Mathematical 
Variable Symbol 

TIME t 

H2 At 

TFIN T 
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a . Scaler Computer Variables, contd 


Computer . 
Variable 

Mathematical 

Symbol 

Description 


CK 

K 

8 

Weighting factor in the performance index. 


FK 

K 

Weighting factor in the performance index. 


EPS 

€ 

Weighting factor in the performance index. 


WEIGT 

k 

Weighting factor in the performance index. 


XMAX 

d 

X 

Translation limit of the cab motion in the x direction. 

YMAX 

d 

Translation limit of the cab motion in the y direction. 


y . 



ZMAX 

d z' ■ ■ 

Translation limit of the cab motion in the z direction. 

f . . . . .. 

2 ; 


* 

FN<t>RM 

g 

max 

Normalization factor used in the approximation of P y . 

RCA 


= COS <t>A . 

. . . 

APHI 

*A 

Aircraft roll angle. 


RSA 


= sin 0 . . . 

A • - i 

' • . 

PCA 


- cos 0 A . 

■ « • ■ 

ATHT 

0 A 

Aircraft pitch angle. 


PSA 


= sin e A . 


YCA 


= COS lf)^. 


APSI 

^A 

Aircraft yaw angle. 


YSA 


=' sin ip 

J 

RNUM 

RDEN 


= sin O ^ cos - cos 0 A sin 0^ sinif)^. ' 

= cos <J> A cos ip^ + sinO ^ sin 0 A sin 


PNUM 


= sin e A . 

> See (2.4) 

*' 

PDEN 


= cos0 A cos i /, A . 


YARG 


= cos 0 A sin^ A . 

i 

IMODE 


Index for control of subroutine computational mode. 
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Computer Mathematical 
Variable Symbol 


Scaler Computer Variables, contd 
Description 


T 

t 

Current time. 


H 

At 

Time increment used for numerical integration. 

IFRST 


Index used to initialize the numerical integration. 

TF2 


= T 2 . 1 


TF3 


= T 3 . 

► See TFIN 

TF4 


J 


CPX 

X 

c 

Cab inertial displacement command in x direction. 

CPY 

y c 

Cab inertial displacement command in y direction. 

CPZ 

X 

c 

Cab inertial displacement command in z direction. 

CVX 

• 

X 

c 

Cab inertial velocity in x direction. 

CVY 

• 

y c 

Cab inertial velocity in y direction. 

CVZ 

• 

z 

c 

Cab inertial velocity in z direction 

CPHI 

0 

c 

Cab roll angle 

command. 

CTHT 

e c 

Cab pitch angle command. 

CPSI 

h 

Cab yaw angle 

• command. 

GAMM1 

*x 

x *• component of the new costate variable y. 

GAMM2 

y 

r y 

y - component of the new costate variable y. 

GAMM3 


z “ component of the new cosfate variable y. 

ATO 


- tan 


APHID 

• 

0 

=P A - (q^ cos 0 A - sin 0 tanfy^; roll gimbal rate. (See 2 .7) 

PA 

pa 

Aircraft roll rate. 

QA 

<A 

Aircraft pitch rate. 

RA 

r A 

Aircraft yaw rate. 

ATHTD 

®A 

- ( q A c «<l>A 

-r^sin0^)/ cosip^; pitch gimbal rate. 

APSID 

$A 

= q A sin0 A + 

cos 0 A '/ yaw gimbal rate. 
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Computer 

Variable 

Mathematical 

Symbol 

Description 

AXFM 

• 'i 

X A 

Aircraft inertial acceleration in x direction. 

AYFM 

• • 

y A 

Aircraft inertial acceleration in y direction. 

AZFM 

• • 

Z A 

Aircraft inertial acceleration in z direction. . , 

CTO 


= tan 0 . 

ABETA2 



CK2 


_ 2 

~ V .. . . . * .. ...r" • 

PV2 


2 

= p . ' 

V 

XXX 



I 


Do loop index. 

CAX 

• • 

X 

c 

x - component of the cab acceleration. 

CAY 

• • 

y c 

y - component of the cab acceleration. 

CAZ 

• • 

z 

X 

z - component of the cab acceleration. A 

CAZZ 

a (3) 
c 

" z c 9 * 

GAM1 


• 

\ 

>■* 

ii 

GAM2 


- v y / k * 

GAM3 


= / k . 

PQR2 

2 

W A 

Magnitude squared of the aircraft angular rate. 

FK2 

K 2 

a) 

Square of the weighting factor in the performance index. . 

GAMMM2 


= (yA> 2 . 

PQRGAM 


= • y / k • 
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Scaler Computer Variables, contd 

Computer 

Variable 

Mathematical 

Symbol 


Description 

OOO 


Ci/ to A * 


PC 

P c 

Cab roll rate. 

QC 

% 

Cab pitch rate. 

RC 

r 

c 

Cab yaw rate. 

CPHID 

• 

<t> 

c 

Derivative of the cab roll angle command. 

CTHID 

• 

0 

C 

Derivative of the cab pitch angle command. 

CPSID 

K 

Derivative of the cab yaw angle command. 



Array Computer Variables 

Computer 

Variable 

Mathematical 

Symbol 

Dimension 

Description 

VC 


12 

Simulator control system state vector. 

YC 


9 

Derivative of the simulator control system state vector, 

AF 

°A 

3 

Specific force acting on the aircraft. 

A ALFA 

“A 

3 

Sensed specific force acting on the aircraft. 

A BETA 

V 

3 

Unbiased sensed specific force acting on the aircraft. 

ACO 


3 

= cos X A - 

ASO 


3 

= sin X A - 

TFM 

c(x A ) 

3x3 

Coordinate transformation of the aircraft Euler angles. 

CCO 


3 

= cos X c « 

CSO 


3 

= sinx c - 

CFM 

c(x c ) 

3x3 

Coordinate transformation of the cab Euler angles. 

CALF A 

— * 
*c 

3 

Sensed specific force acting on the cab. 
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Array Computer Variables, contd 


Computer 

Variable 

Mathematical 

Symbol 

Dimension 

Description 

PV 


3 

Costate of the velocity error. 

ROA 


3 

= c (> C )S A - 

ROB 


3 

= c '(x c )v 

W 

-4 

a 

3 

Error in the specific force. 

ETA 

V 

3 

Intermediate variable used in the angular washout 
computation. 

D 


3 

Diagonal elements of the matrix D in the penalty 

function L (7 ). 

c 


14? 



o o o o o 


Listing of Subroutine WASHFL 


SUBROUTINE WASHFL ( I MODE , H, APH I , ATHT , APS I , AXFM , A YFM , AZFM , 
* PA ,QA »R A , CPX »CPY ,CPZ » CPH I , CTHT ,CPS I > 

DIMENSION VC <12 ) ,YC(9) 

DIMENSION AF < 3 ) ,AALFA(3) , ABETA ( 3 ) , ACO( 3 ) , ASO ( 3 ) ♦ tFM ( 3 , 3 ) 
DIMENSION CCO(3) ,CS0(3) ,CFM( 3,3) ,CALFA(3) 

DIMENSION PV ( 3 ) ,R0A(3) , ROB ( 3 J , W < 3 ) , ET A ( 3 ) , D ( 3 ) 


DATA TFIN/5.0 / 
DATA CK/C.O / 
DATA FK/C.O / 
DATA EPS/C. 16 / 
DATA WEIGT/1.0 / 
DATA XMAX/10.0 / 
DATA YMAX/10.0 / 
DATA ZMAX/10.0 / 


DATA FNORM/3 50 • 0 / 

DOUBLE PRECISION TIME.H2 
RCA= COS( APHI ) 

RSA= SINIAPHI) 

PC A= COS I ATHT ) 

PSA= S I N 1 ATHT ) 

YCA= C0S( APSI ) 

Y S A = S IN ( APS I ) 

RNUM = — RCA*PSA*YSA + RSA*YCA 
RDEN = RSA*PSA*YSA + RCA*YCA 
PNUM = PSA 
PDEN » PC A* YC A 
YARG = PC A*Y SA 
APHI = ATAN2(RNUM*RDEN) 

ATHT = ATAN2IPNUM.PDEN) 

APSI = ASIN(YARG) 

IF (IMODE) 100C, 3000, 2000 


MAYBE 

MAYBE 

MAYBE 

MAYBE 

MAYBE 

MAYBE 

MAYBE 

MAYBE 

MAYBE 

MAYBE 

MAYBE 

MAYBE 

MAYBE 

MAYBE 


***#*##$*$$**# 

INITIALIZATION 

#*#*#$**★#$*** 


1000 CONTINUE 
T = 0.0 

TIME = DBLE(T) 

T = SNGL ( TIME) 

H2 = DBLE(H) 

IFRST = 0 

TF2 = TFIN**2 
TF3 = TFIN**3 
TF4 = TF2**2 
D( 1) = 1.0/(XMAX**2) 
0(2) = 1.0/(YMAX**2> 
D ( 3 ) = 1 .0/ ( ZMAX**2 ) 
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Listing of Subroutine WASHFL (Continued) 


CPX = 

0.0 

CPY- = 

C.O 

CPZ = 

0.0 

cvx = 

0.0 

CVY = 

0.0 

cvz = 

0.0 

CPH I 

= APHI 

CTHT 

= ATHT 

CPSI 

= APSI 

GAMM1 

= 0.0 

GAMM2 

= 0.0 

GAMM3 

= 0.0 

GO TO 

3000 


DEFINING AIRCRAFT VARIABLES 


2000 CONTINUE 

ACO(l) = COS(APHI) 

ACO I 2 ) = COS(ATHT) 

ACO ( 3 ) = COS ( APS I ) 

A SO ( I ) = SIN< APHI ) 

A SO { 2 ) = SIN(ATHT) 

A SO ( 3 ) = SINIAPSI ) 

ATO = TAN(APSI) 

APHID = PA - (CA*AC0<1) - RA* ASO I 1 ) ) *ATO 
ATHTD = <QA*AC0(1) - RA*AS0 (1) ) / ACO ( 3 ) 
APSID = QA* ASC ( 1 ) + R A* ACO ( 1 ) 


C 


c 


TFM ( 1 , 1 ) 
TFMI2.1) 
TFM I 3 , 1 ) 
TFM( 1,2) 
TFM ( 2 * 2 ) 
TFM ( 3 » 2 ) 
TFM I 1,3) 
TFM ( 2 , 3 ) 
TFM (3, 3) 


ACO ( 2 ) * ACO ( 3 ) 

A SO ( 1 ) ♦ A SO t 2 ) - ACO( 1 )*AS0(3)*AC0(2) 
ACO ( 1 ) * A SO ( 2 ) + ASOI 1 ) *ASO I 3 ) *ACO < 2 ) 
A SO ( 3 ) 

ACCI 1 ) * ACO ( 3 ) 

-A SO ( 1 ) * ACO ( 3 ) 

-ACOI3)*ASO(2) 

ASOI 1)*AC0<2) + ACCI 1)*AS0(3)*AS0(2) 
ACO I 1 ) ♦ACO I 2 ) - ASO(l)*ASO(3)*ASOI2> 


AF(l) = AXFM 
AFI2) = A YFM 
AF I 3 ) = AZFM. - 32.2 

AALFA(l) = TFM ( 1 , 1 ) $ AF ( 1 ) + TFM I 1,2) *AF I 2 ) + TFMI 1 

AALFAI2) = TFM I 2 , 1 ) *AF I 1 ) + TFM I 2 , 2 ) * AF I 2 ) + TFM I 2 

AALFAI3) = TFM ( 3 , 1 ) * AF 1 1 ) + TFM I 3 , 2 ) *AF ( 2 ) + TFM I 3 


3 ) * AF I 3 ) 
3 ) * AF I 3 ) 
3 )*AF( 3) 
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Listing of Subroutine WASHFL (Continued) 


ABETA ( 1) .= AALFA(l) ; 

ABETA ( 2 ) = A ALFA ( 2 ) 

A B E T A ( 3 ) = AALF A ( 3 ) ♦ 32.2 




DEFINING CAB VARIABLES 
$$**#$$*$*#*#**«#«$**$ 



CCO( 1 ) = COS ( CPH I ) 
CC0I2) = COS(CTHT) 
CCO ( 3 ) = COSICPSI ) 
CSOI 1 ) = SINICPHI ) 
CSO ( 2 ) = SIN(CTHT) 
C SO ( 3 ) = SIN (CPS I) 
CTOff = TANICPSI ) 


CFM( 1,1) 
CFM ( 2 i 1 ) 
CFM ( 3 » 1 ) 
CFM ( 1 » 2 ) 
CFM ( 2 f 2 ) 
CFM ( 3 » 2 I 
CFM ( 1 » 3 ) 
CFM ( 2 » 3 ) 
CFM ( 3 i 3 ) 


CCO ( 2 ) *CCO ( 3 ) 

CSOI l)*CSO(2) - CCO(l)*CSO(3)*CCO(2) 
CCO ( 1 ) *CSO ( 2 ) + CS0T1 >*CSQ(3)*CCO(2) 
CSC (3) 

CCC( 1 )*CCO(3) 

-CSO( l)*CCO(3> 

-CCO ( 3 ) *C SO ( 2 ) 

CSO(l)*CCO(2) + CCOm*CSO(3)*CSOl2> 
CCO ( I ) *CCu ( 2 ) - CSO(1)*CSO(3)*CSO(2) 


COMPUTES TRANSLATIONAL WASHOUT CONTROLS ; , c . , 

**#*$**$$****#* * **** * **** * ***** * $ 

1 • ’ ‘ } 

PV(1) = -EPS*D(I )*(CPX*TF2/2.0 +CVX*TF3/3.0 +AXFM*TF4/(8.0*FN0RM) ) 
PV ( 2 ) = -EPS*D(2)*(CPY*TF2/2.0 +CVY*TF3/3.0 +AYFM*TFA/(8.0*FNORM) ) 
P V ( 3 ) = -EPS*D(3)*(CPZ*TF2/2.0 +CVZ*TF3/3.0 .+ AZFM*TF4/ ( 8.0*FN0RM ) ) 
ROA(l) = AALFA(l)*CFMll,l) + AALF A ( 2 ) *CFM ( 2 » 1 ) + AALF A ( 3 ) *CFM ( 3 , 1 ) 
RO A ( 2 ) = AALFAIl )*CFMU,2) + AALFA ( 2 ) *CFM ( 2*2) + AALF A ( 3 ) *CFM ( 3 , 2 ) 

RO A ( 3 ) = AALFAm*CFMU,3) + A ALF A.( 2 ) *CFM ( 2 , 3 ) + AALF A I 3 ) *CFM l 3 , 3 ) 

ABETA2 = ABET A ( 1 ) **2 + ABETA(2)**2 + ABETA(3)**2 . 

C 

IF (CK.EQ.O.O) GO TO 100 
CK2 = CK**2 

PV2 = PV ( 1 ) **2 + PV ( 2 ) **2 + PV<3)**2 

ROB ( 1 ) = ABETA ( 1 ) *CFM ( 1 » 1 ) + ABET A ( 2 ) *CFM ( 2 „1 ) +, ABET A ( 3 ) *CFM ( 3 , 1 ) 

ROB ( 2 ) = ABETA(1)*CFM(1,2) + ABETA < 2 ) *CFM ( 2 ,2 ) + ABET A ( 3 ) *CFM ( 3 , 2 > 

ROB ( 3 ) = ABETA(1)*CFM(1,3) + ABET A ( 2 ) *CFM ( 2 , 3 ) + ABET A ( 3 ) *CFM ( 3 , 3 ) 

XXX = ABETA2*PV2 + CK2 + 2.*CK * ( PV (1) *ROB ( 1 ) +PV ( 2 ) *ROB ( 2 ) 

* ;+PV(3)*R0B(3) ) 

XXX = SORT ( XXX ) 
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Listing of Subroutine WASHFL (Continued) 


C 


DO 5 1=1,3 

ETA ( I ) = - ( ( 1 . C + CK/(-CK+XXX) ) ALFA ( I ) + CK*ABETA( I ) ) / ABET A2 
5 W(I) = (1.0 - CK/XXX)*(ABETA2*PV( I ) ♦ CK*R0B(I)> + ROA(I) - AF ( I ) 
GO TO 200 
100 CONTINUE 


DC 6 1=1,3 

ETA ( I ) = -AALFA ( I ) /ABETA2 
6 W(I) = ABET A2*PV ( I ) + ROA(I) - AF(I) 


200 CONTINUE 
C AX = AXFM 
CAY = AYFM 
CAZ = A ZFM 
CAZZ = CAZ 
C ALFA ( 1 ) = 
CALFA ( 2 ) = 
C ALFA ( 3 ) = 


+ W(l) 

+ W(2> 

♦ W<3) 

- 3 2.2 

CFM ( 1 » 1 ) *CAX 
CFM ( 2,1) *CAX 
CFM ( 3 , 1 ) *CAX 


+ CFM ( 1 , 2 ) *C AY 
+ CFM ( 2 , 2 ) *CAY 
+ CFM { 3 , 2 ) *C AY 


♦ CFM( 1 » 3 ) *C AZ Z 

♦ CFM (2,3) #CAZZ 

♦ CFM ( 3 , 3 ) *C AZZ 


«**$#$«*###$*#*$«#$$$*#****#**♦#* 


COMPUTES ANGULAR WASHOUT CONTROLS 
$*$*$*********♦**************$**$ 


G AM 1 = GAMM1/WEIGT 
GAM2 = GAMM2/WE IGT 
GAM3 = GAMM3/WE IGT 
PGR2 = P A*PA + QA*QA ♦ RA*RA 
C 

IF (FK.EQ.O.O) GO TO 700 ! 

FK2 = FK*FK 

GAMMM2 = GAM1**2 ♦ GAM2**2 ♦ GAM3**2 
PQRGAM = PA*GAM1 ♦ Q A*GAM2 + RA*G AM3 
000 = PQR2*GAMMM2 + 2 . 0*FK*PCRGAM + FK2 
000 = SORT (000) 

PC -PA + (1.0 - FK/Q00)*(PQR2*GAMl + FK*PA) 
QC = QA + (1.0 - FK/G00)*(PQR2*GAM2 + FK*QA> 

RC = RA + (1.0 - FK/QOO)*(POR2*GAM3 + FK*RA ) 

GO TO 800 
700 CONTINUE 
C 

PC = PA + PQR2*GAM1 

GC = QA + PQR2*GAM2 

RC = RA + PQR2*GAM3 

800 CONTINUE 
C 

CPHID =PC - QC*CCO( 1 )*CTO + RC*CS0(1)*CT0 
CTHTD =QC*CCO( 1 ) /CCO ( 3 ) - RC*C SO ( 1 > /CCO ( 3 ) 
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Listing of Subroutine WASHFL (Continued) 


CPSIO =QC*C SO ( 1 ) ♦ RC*CCQ ( 1 ) 


**#*##*$*$*#*$*$**$***#*$#«$##$$## 
COMPUTES INPUTS TO THE INTEGRATORS 


YC ( I ) 

= 

C AX 

YC (2 ) 

= 

CAY 

YC ( 3 ) 

= 

CAZ 

YC ( 4 ) 

= 

CPHID 

YC l 5 ) 

= 

CTHTD 

YC C 6 ) 

= 

CPSIO 

YC ( 7 ) 

SB 

CALFA(3)*ETA!2) 

YC ( 8 } 

= 

CALFA! 1 >*ETA!3) 

YC ( 9 ) 

= 

CALFA ( 2 ) *ETA { 1 ) 

VC(1) 

— 

CPX 

VC (2) 

= 

CVX 

VC ( 3 ) 

= 

CPY 

VC (4) 

3= 

CVY 

VC ( 5 ) 

= 

CPZ 

VC ( 6 ) 

= 

CVZ 

VC ( 7 ) 

SB 

C PH I 

VC (8) 

= 

CTHT 

VC ( 9 ) 

= 

CPSI 

V C ( 10 > 


= GAMM1 

VC C 1 1 ) 


= GAMM2 

VC( 12) 


= GAMM3 


CALFA ( 2 ) *ETA ( 3 ) 
CALFA ! 3 ) *ETA ( 1 ) 
CALFA { 1 ) #ETA ( 2 ) 


PERFORMING INTEGRATION 


C 


CALL 

GINTRtT, H, 

CPX = 

VC( 1) 

CPY = 

VC( 3) 

CPZ = 

VC! 5) 

CVX = 

VC (2) 

CVY = 

VC (4) 

CVZ = 

VC ( 6 ) 

CPH I 

= VC (7) 

CTHT 

= VC (8) 

CPSI 

= VC ( 9 ) 

GAMM1 

= VC! 10) 

GAMM2 

= VC! 11) 


YC, IFRST) 
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Listing of Subroutine WASHFL (Continued) 


GAMM3 = VC ( 12) . ' 

c \ , .. 

TIME = TIME + H2 
T = SNGL (TIME) 

C 

3000 CONTINUE 
RETURN 
END 



Listing of Subroutine GINTR 


SUBROUTINE GINTR ( T f H, Y, YC, IFRST) 
DIMENSION Y(12),YD(12),YC!9t,XXX(12> 
DIMENSION ST ATE ( 1 2 ) , STATET ( 12 ) 
COMMUN/DOUBLE/ STATE, STATET 
DOUBLE PRECISION STATE, STATET, DTT , DBLT 
T MAX = T + H 

IF ( I FRST.NE . 0 ) GO TO 900 
IFRST = 1 
EPSIL = l.E-4 
CALL DRPRCIYQ, Y, YC i 
900 CONTINUE 
101 CONTINUE 

DBLT = DBLE(T) 

DTT = DBLT ♦ 0.5D0 * DBLE(H) 

TT = SNGL ( DTT ) 

DO 105 1=1,12 

STATE ( I ) = DBLE ( Y ( I ) > 

XXX( I ) = 0.5*H*YD< I ) 

STATE T(I) = STATE ( I t+DBLE (XXXC I ) ) 

Y ( I )=SNGL( STATET! I ) > 

105 CONTINUE 

CALL ORPRC ( YD , Y, YC ) 

DO 205 1=1,12 
XXX! I ) = H* YO ( I ) 

STATE { I > = STATE ( I )+DBLE(XXX( I ) ) 

Y ( I ) = SNGL ( STATE (I) ) 

205 CONTINUE 

DBLT = DBLT + DBLE(H) I 

T = SNGL ( DBLT ) 

CALL DRPRCIYD, Y, YC ) 

IF ( T. LT . ( TMA X - EPSIL )) GO TO 101 

RETURN 

END 

» 
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Listing of Subroutine DR PRC 


C 

C 


C 

C 


SUBROUTINE ORPRC ( DXP , ZP * VP ) 
DIMENSION DXP C 1 2 ) ,ZP(12> ,VP(9> 


DXP ( 1 ) 

= 

ZP ( 2 ) 

DXP ( 2 ) 

= 

VP( 1 > 

D XP ( 3 ) 

= 

ZP ( A ) 

DXP ( A ) 

- 

VP ( 2 ) 

DXP ( 5 ) 

= 

ZP ( 6 ) 

DXP ( 6 ) 


VP( 3) 

DXP ( 7 ) 

= 

VP ( A ) 

DXP ( 8 ) 

= 

VP ( 5 ) 

DXP (9) 

= 

VP ( 6 ) 


DXP ( 10) = VP ( 7 ) 
DXP (11) = VP ( 8) 
DXP ( 1?) = VP ( 9 ) 


RETURN 

END 


*U.S GOVERNMENT PRINTING OFFICE: 19 73-739-153/49 
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